Data analysis from measurements (PRBS injected to a system -> how to plot the bode of the system response)
3 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Yannick CAIRE
il 20 Mag 2022
Commentato: Yannick CAIRE
il 16 Ott 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
Risposta accettata
Mathieu NOE
il 20 Mag 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 Commenti
Mathieu NOE
il 23 Mag 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
Più risposte (1)
Rajiv Singh
il 31 Mag 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.
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!