Filter löschen
Filter löschen

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

6 Ansichten (letzte 30 Tage)
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

Akzeptierte Antwort

Mathieu NOE
Mathieu NOE am 20 Mai 2022
hello
see my suggestion below (data in attachement)
clc
clearvars
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
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
  3 Kommentare
Mathieu NOE
Mathieu NOE am 23 Mai 2022
hello again
if you need a good quality FRF down to 0.1 Hz resolution , you need NFFT = 10 x Fs at least.
I don't know your Fs but the amont of samples must be larger than NFFT itself... so can be pretty large.
also plot your result with x axis in log scale to compare with the theoretical Bode plot.
all the best
Yannick CAIRE
Yannick CAIRE am 24 Mai 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.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (1)

Rajiv Singh
Rajiv Singh am 31 Mai 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 Kommentar
Yannick CAIRE
Yannick CAIRE am 16 Okt. 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.

Melden Sie sich an, um zu kommentieren.

Produkte


Version

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by