# Data analysis from measurements (PRBS injected to a system -> how to plot the bode of the system response)

4 views (last 30 days)
Yannick CAIRE on 20 May 2022
Commented: Yannick CAIRE on 16 Oct 2022
Hello,
I have data (measurements from a synchronous machine), a PRBS has been injected to the excitation system (AVR setpoint) and the output measurements have been collected (Active Power for example). So I have the perturbed setpoint (Set point of the regulator with the PRBS added to this SP, which is my input u), and I have the ouput y (Active Power P).
Now I need to trace the bode of P/SP (y/u). I'm mostly interseted in the low frequencies from 0.05 Hz to 20Hz (max 100Hz).
I would like to know, how I can proceed in Matlab, do I need a particular Toolbox ?
Any advice would be highly appreciated.
Thanks a lot
##### 2 CommentsShowHide 1 older comment
Yannick CAIRE on 20 May 2022
Thanks a lot, do you have an example similar to my case ? Thanks

Mathieu NOE on 20 May 2022
hello
see my suggestion below (data in attachement)
clc
clearvars
x = transpose(data.x); %input
y = transpose(data.y); %output
fs = data.fs; % sampling frequency
NFFT = 2048;
NOVERLAP = round(0.75*NFFT); % 75 % overlap
%% solution 1 with tfestimate (requires Signal Processing Tbx)
% [Txy,F] = tfestimate(x,y,hanning(NFFT),NOVERLAP,NFFT,fs);
%% alternative with supplied sub function
[Txy,Cxy,F] = mytfe_and_coh(x,y,NFFT,fs,hanning(NFFT),NOVERLAP); % t = transfer function (complex), Cxy = coherence, F = freq vector
% bode plots
figure(1),
subplot(2,1,1),plot(F,20*log10(abs(Txy)));grid
subplot(2,1,2),plot(F,180/pi*(angle(Txy)));grid
%%%%%%%%%%%%%%%%%%%%%%%
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
Yannick CAIRE on 24 May 2022
Thanks.
Yes fs is 200Hz (0.005s period) and I have a large amount of points (24000).
It's gettting better with a NFFT around 6000 for example, but I think I can still improve it.

Rajiv Singh on 31 May 2022
This could be seen as an exercise in system identification. Use functions like TFEST,SSEST, ARX to fit a model to the (u,y) data. Then call BODE or FREQRESP on the model to get the frequency response.
Requires System Identification Toolbox.
Yannick CAIRE on 16 Oct 2022
Thanks, do you have an example using these functions ? Then I could start from somewhere, I have the input and ouput data captured at 5ms sampling interval, using a PRBS. Input is the PRBS + the setpoint of the system.

### Categories

Find more on Spectral Analysis in Help Center and File Exchange

R2021b

### Community Treasure Hunt

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

Start Hunting!