Determine of coherence using Transfer function (Bode plot)

19 views (last 30 days)
Hello everyone,
I was trying to understand gain and phase margins for the flow fluctuations in my domain. During the process, I came across the following blog (https://in.mathworks.com/matlabcentral/answers/1628890-determine-the-bode-diagram-or-a-transfer-function-with-input-output-data-without-tfest ), where I could able to understand the gain and phase margins + coherence function plot based on the transfer function.
I got the code from the resepective link. I calculated the transfer function based on power spectra for my input and output data and plotted gain and phase margin as shown in image below.
However, coherence is maintained one, without fluctuations, that is puzzling me -- is that mean, energy is always higher in my signal ? Can someone share some insights and correct me if I am wrong.
Thank you
clc; clear all;
clc
clearvars
%importdata('velocityW.csv');
load('velocityW.mat');
input = ans(:,2);
output = ans(:,3);
time = ans(:,1); %*10^(-6);
Fs = 1/mean(diff(time)); % Sampling Frequency
nfft = 10000;
window = hanning(nfft);
noverlap = round(0.75*nfft);
[Txy,Cxy,Fv] = tfe_and_coh(input,output,nfft,Fs,window,noverlap);
figure
subplot(3,1,1)
semilogx(Fv, 20*log10(abs(Txy)))
title('Amplitude')
ylabel('dB')
subplot(3,1,2)
semilogx(Fv, angle(Txy)*180/pi)
title('Phase')
ylabel('°')
subplot(3,1,3)
semilogx(Fv, Cxy)
title('Coherence')
%%%%%%%%%%%%%%%%%%%%%%%
function [Txy,Cxy,f] = tfe_and_coh(x,y,nfft,Fs,window,noverlap)
% Transfer Function and Coherence Estimate
% compute PSD and CSD
window = window(:);
n = length(x); % Number of data points
nwind = length(window); % length of window
if n < nwind % zero-pad x , y if length is less than the window length
x(nwind)=0;
y(nwind)=0;
n=nwind;
end
x = x(:); % Make sure x is a column vector
y = y(:); % Make sure y is a column vector
k = fix((n-noverlap)/(nwind-noverlap)); % Number of windows
% (k = fix(n/nwind) for noverlap=0)
index = 1:nwind;
Pxx = zeros(nfft,1);
Pyy = zeros(nfft,1);
Pxy = zeros(nfft,1);
for i=1:k
xw = window.*x(index);
yw = window.*y(index);
index = index + (nwind - noverlap);
Xx = fft(xw,nfft);
Yy = fft(yw,nfft);
Xx2 = abs(Xx).^2;
Yy2 = abs(Yy).^2;
Xy2 = Yy.*conj(Xx);
Pxx = Pxx + Xx2;
Pyy = Pyy + Yy2;
Pxy = Pxy + Xy2;
end
% Select first half
if ~any(any(imag([x y])~=0)) % if x and y are not complex
if rem(nfft,2) % nfft odd
select = [1:(nfft+1)/2];
else
select = [1:nfft/2+1]; % include DC AND Nyquist
end
Pxx = Pxx(select);
Pyy = Pyy(select);
Pxy = Pxy(select);
else
select = 1:nfft;
end
Txy = Pxy ./ Pxx; % transfer function estimate
Cxy = (abs(Pxy).^2)./(Pxx.*Pyy); % coherence function estimate
f = (select - 1)'*Fs/nfft;
end

Answers (1)

William Rose
William Rose on 16 Dec 2024
The coherence equals one in your example because you estimated the spectra Pxx, Pyy, and Pxy using a single segment. If you use multiple segments, the coherence esitmate will have be between 0 and 1. I will explain this more below, but first:
The time vector in your data is sampled at a non-uniform rate. See plot of delta t versus sample number below.
The plot shows that there are handful of points sampled fast, then about 60% of the points are sampled relatively slowly, then the last 40% (approximately) of points are sampled about 17 times faster (dt is 17 times smaller) than the first 60%. Therefore you should resample to a uniform rate, with interp1(), before estimating spectra. Or maybe the time vector in your data file is wrong.
Let's pretend the data is sampled at a uniform rate, with dt=1, i.e. fs=1/dt=1.
Estimate coherence the way you did, with a single segment that includes all the data.
load('velocityW');
x=ans(:,2); y=ans(:,3);
N=length(x); dt=1;
Pxx1=cpsd(x,x,N,0,N); % estimate Pxx with a single full-length segment
Pyy1=cpsd(y,y,N,0,N); % estimate Pyy with a single full-length segment
Pxy1=cpsd(x,y,N,0,N); % estimate Pxy with a single full-length segment
Coh1=abs(Pxy1).^2./(Pxx1.*Pyy1); % coherence estimate with a single segment
f1=(0:(N-1)/2)/(dt*N); % frequency vector
Next: estimate coherence with segments that are one-eighth of full length (or just a bit shorter than 1/8, so that 8 segments will fit), plus 7 segments of equal length that overlap the first 8 by half on each side. Thus there will be 15 total segments, each with length <= 1/8 of the full length sequence. This approach uses all or almost all the data, and reduces the variance in the power spectral estimate, compared to the estimate obtained with a single full-length segment.
Nseg=floor(N/8); % one-eighth length, or a bit less
if mod(Nseg,2)==1,Nseg=Nseg-1; end % if Nseg odd, decrement by 1 to make Nseg even
Pxx8=cpsd(x,x,Nseg,Nseg/2,Nseg); % estimate Pxx with 15 half-overlapped segments
Pyy8=cpsd(y,y,Nseg,Nseg/2,Nseg); % estimate Pyy with 15 segments
Pxy8=cpsd(x,y,Nseg,Nseg/2,Nseg); % estimate Pxy with 15 segments
Coh8=abs(Pxy8).^2./(Pxx8.*Pyy8); % coherence estimate
f8=(0:Nseg/2)/(dt*Nseg); % frequency vector
Plot results
figure
plot(f1,Coh1,'-r.',f8,Coh8,'-b.')
grid on; xlabel('Frequency'); ylabel('\gamma^2'); title('Coherence Estimates')
legend('1 segment','15 segments');
The code above produces the figure below:
The figure shows that when one estimates coherence with one segment, the coherence equals 1 at all frequencies. But if one uses multiple segments, the estimate of coherence is more realistic, and is between 0 and 1. Frequency goes from 0 to 0.5 in the plot above, because we assumed sampling rate fs=1. Thefore the coherence spectra go up to the Nyquist frequency.
Computing coherence with a single segment is kind of like like estimating Pearson correlation (r-squared) for the best-fit straight line, for a data set with just two points. The line fits the data perfectly when there are only two points, so r^2=1.000. Not a perfect analogy, but close.
  2 Comments
Kumaresh Kumaresh
Kumaresh Kumaresh on 16 Dec 2024
Thank you for clarifying here. Yes you are right. During my simulation, first half was calculated very slowly (good sinusoidal behavior in flow pattern), whereas second half was calculated with faster time step. The reason why you observe disorderness and strange behavior in the second half profile. So, if I could able to extract the time data in a uniform manner (clear signal), I believe my data will be reasonable to interpret.
Thank you once again. I will update my results here in some days (after extracting the clean signal) to discuss with you once again.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!