Determine of coherence using Transfer function (Bode plot)

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

Answer 1

[William Rose]

@Kumaresh Kumaresh,

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.

[Kumaresh Kumaresh]

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.

[William Rose]

@Kumaresh Kumaresh, you're welcome. I look forward to seeing your results.

[Kumaresh Kumaresh]

Hello Mr. William Rose,

Hope you are doing good. Based on our previous conversations, I would like to share with you my plot with some questions.

here is the Bode plot for pressure variable

For velocity

Here is my final graph explained for pressure and velocity.

Pressure is not the dominating parameter is my study, as we could see at higher frequencies, amplitude and phase becomes zero, whereas the velocity dominates at higher frequency, for input frequencies of 100 and 200Hz, when the amplitude becomes zero, the phase shift is higher oscillating between +180 to -180, proving unstability.

Kindly correct me if I am wrong and tell me if I missed something in the plot.

In addition, I would like to plot Nyquist or Nichols plot. Because in Nyquist or Nichols plots, the information will be in single chart + easier to determine stability, however frequency is not explicitly called out.

And I read that Nyquist plot easily handles time delay, which I am curious to plot and reading more about it.

Sorry for too many questions.

Share your suggestions please.

Thank you

[William Rose]

@Kumaresh Kumaresh,

Your latest post shows plots of magnitude, phase and coherence for pressure and for velocity. Bode plots and coherence are defined between two variables. For the "pressure" plots, what is in the numerator and what is in the denominator? The same question applies to the "velocity" plots.

You say "the phase shift is higher oscillating between +180 to -180, proving unstability". That is true if this is a plot of dimensionless gain in the open loop condition. The gain and phase margin rules for stability are applicable for plots of open loop gain when the input and output have the same units. Did you open the feedback loop in this experiment?

Nichols plots and Nyquist plots are also for open loop data.

An example of opening the feedback loop is in the study of carotid baroreceptor control of arterial blood pressure. You do surgery to separate the blood pressure receptor region from the rest of the circulation, while preserving the nerves carrying pressure information to the brain. This opens the feedback loop. I have done those experiments in multiple species, and so have many others. For example. see figure below, from Prof. Kenji Sunagawa et al., 2019. The top pressure in the figure is the input (carotid sinus pressure, CSP), the pressure below it is the output (arterial pressure, AP). The open loop gain is always less than 1, so the system is stable when the loop is closed. And in the intact animal, the loop IS closed.

If you provide the data and code used to make the figures, and explain if the data is open loop or not, people will be able to offer more useful advice.

[Kumaresh Kumaresh]

Hello Mr. William Rose,

Thank you for your response. My work is related to computational simulation of aircraft gas turbine combustor. Based on the input conditions (in the form of sine wave) applied, output distorted sine wave is obtained. Based on the physics, pressure and velocity are two critical parameters, which are examined.

So as per your question, I didn't apply numerator or denominator specifically like other traditional problems. No feedback loop analysis done for this open-loop data. Thank you for explaining with an example. There might be no possibility to compute a valid transfer function in a given frequency range (as no information from the input is observed in the outputs). With the given problem I have, I approached Bode Plot to understand the behavior in the frequency domain and examine the stability.

Here is my code (Gain_PhaseMargin.m) and data (test.mat) attached in this thread.

In test.mat data file,

1 column - time, 2 column - input, all other columns - output, can be compared with input singal separately.

Kindly share your valuable insights.

Thank you

[William Rose]

@Kumaresh Kumaresh,

I downloaded your code and your data. Script Gain_PhaseMargin.m runs without error.

After the code ran, I computed the sampling interval between each sample and plotted the results.

load('test.mat');
time=ans(:,1);
dt=diff(time); N=length(time);     % compute sampling interval
disp([min(dt),mean(dt),max(dt)])  % display minimum, mean, maximum

   1.0e-05 *

    0.0456    0.4007    0.5000

figure;
subplot(211), plot(1:N-1,dt,'.k') 
grid on; title('Sampling Interval vs. Sample Number')
xlabel('Sample Number'); ylabel('\Delta T (s)')
subplot(212), plot(time(2:end),dt,'.k')
grid on; title('Sampling Interval vs. Time')
xlabel('Time (s)'); ylabel('\Delta T (s)')

The sampling interval varies from 0.45 microseconds to 5 microseconds. The frequencies of the discrete Fourier transform are computed with the assumption of a constant sampling rate. The sampling rate here is not constant: it varies by more than a factor of ten. Therefore the frequency analysis is not reliable. You may be able to fix this by intepolationg to the fastest rate. The results of that interpolation will only be valid if the signal does not have power at frequencies faster than one half of the slowest sampling frequency in the data.

I assume that the data columns which you call "input" and "output" are pressure and velocity. I do not know which is which. In any case, it is impossible to understand anything about stability of this system from this data. Gain and phase margin analysis is for when you have dimensionlees transfer funcitons in the open loop conditions. Your transfer functions are not dimensionless (they have units of pressure/velocity, or the inverse) and you said the data is not in the open loop condition.

You have 5000 data points. You set nfft=4000 and noverlap=1000. Therefore you will have two windows worth of data when you compute the spectra, and the windows overlap for 3/4 of their length. Data from windows that overlap so highly have a high correlation, and therefore they do not reduce the variance of the spectral estimates as much as windows that do not overlap, or which overlap by a smaller fraction of their length. Therefore, in order to get spectral esitmates that less variability, I recommend using a smaller value for nfft, and noverlap=nfft/2. The drawback of this is that the lowest frequency in the spectra will be higher, and the frequency resolution will be lower.

[Kumaresh Kumaresh]

Sorry for replying late. Thank you for your comments.

Transfer function for pressure = p'_combustor / p'_inlet; where p' is pressure fluctuation.

Transfer function for velocity = u'_combustor / u'_inlet; where v' is velocity fluctuation.

So my gain is dimensionless.

[William Rose]

@Kumaresh Kumaresh,

Do you mean that the two variables you provided in test.mat are already dimensionless transfer functions? If so, how did you calculate them? If so, then we have been calculating the ratio of two transfer functions, which one coudl think of as the trasfer funciton from one to another. That is a concept that I have never come across before.

[Kumaresh Kumaresh]

@William Rose

p'_combustor and p'_inlet has the unit of Pascal. When they are implemented in transfer function, the transfer function becomes dimentionsless.

Transfer function for pressure = p'_combustor / p'_inlet = Pa / Pa

p'_inlet is the one of inputs in the form of sine wave in my computational simulation and p'_combustor is the calculated output.

Thank you