hello
see example below
from the Bode plot you can easily extract the plant dc gain (here 10) and cut off frequency (fc = 50 Hz here, corresponding to a -3 dB gain drop vs dc gain and a phase rotation of -45 degrees)
fs = 1000; % sampling frequency (Hz)
dt = 1/fs;
%% generate the test signal
samples = 10000;
% x = rand(samples,1); % random
x = 0.5*(1+sign(randn(samples,1)));% "home made" simplified PRBS
%% plant = 1st order filter
N = 1; % filter order
fc = 50; % Hz cut off frequency
DC_gain = 10; % 1st order plant dc gain
%% generate the filtered output signal
[B,A] = butter(N,2*fc/fs);
B = B*DC_gain;
y = filter(B,A,x);
t = (0:samples-1)*dt;
figure(1)
subplot(2,1,1),plot(t(1:100),x(1:100),'b')
subplot(2,1,2),plot(t(1:100),y(1:100),'r')
%% solution 1 with tfestimate (requires Signal Processing Tbx)
% [Txy,F] = tfestimate(x,y,hanning(NFFT),NOVERLAP,NFFT,fs);
%% alternative with supplied sub function
NFFT = 512;
NOVERLAP = round(0.5*NFFT); % 0 percent overlap
[Txy,Cxy,F] = mytfe_and_coh(x,y,NFFT,fs,ones(NFFT,1),NOVERLAP);
% Txy = transfer function (complex), Cxy = coherence, F = freq vector
% Bode plots
fmax = fs/2.56; % plot only values of freq between 0 and fs/2.56 to avoid drop at Nyquist frequency
figure(2),
subplot(3,1,1),semilogx(F,20*log10(abs(Txy)));
xlim([0 fmax]);
ylabel('Mag (dB)');
subplot(3,1,2),semilogx(F,180/pi*(angle(Txy)));
xlim([0 fmax]);
ylabel('Phase (°)');
subplot(3,1,3),semilogx(F,Cxy);
xlim([0 fmax]);
ylim([0 1]);
xlabel('Frequency (Hz)');
ylabel('Coherence');
%%%%%%%%%%%%%%%%%%%%%%
function [Txy,Cxy,f] = mytfe_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
