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

6 views (last 30 days)
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

Accepted Answer

Mathieu NOE
Mathieu NOE on 20 May 2022
see my suggestion below (data in attachement)
data = load('beam_experiment.mat');
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
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 = 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;
% 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];
select = [1:nfft/2+1]; % include DC AND Nyquist
Pxx = Pxx(select);
Pyy = Pyy(select);
Pxy = Pxy(select);
select = 1:nfft;
Txy = Pxy ./ Pxx; % transfer function estimate
Cxy = (abs(Pxy).^2)./(Pxx.*Pyy); % coherence function estimate
f = (select - 1)'*Fs/nfft;
Yannick CAIRE
Yannick CAIRE on 24 May 2022
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.

Sign in to comment.

More Answers (1)

Rajiv Singh
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.
  1 Comment
Yannick CAIRE
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.

Sign in to comment.




Community Treasure Hunt

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

Start Hunting!

Translated by