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

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

[Star Strider]

System Identification Toolbox

[Yannick CAIRE]

Thanks a lot, do you have an example similar to my case ? Thanks

Answer 1

[Mathieu NOE]

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

[Yannick CAIRE]

Thanks a lot, I tried first with the Method 1 (using Signal Porcessing Toolbox). The result is not far from what I expect, but I need something more precise on the low frequency (0.1Hz to 10Hz is the area of interset).

I probably need to adjust NFFT or another parameter ?

See the attached pictures, Ftestimate is what I get with Method 1, and Bode is the simulation result. So I expect the Signal Processing Toolbox to give me a "similar result", I can see it's almost there but probably something needs to be adjusted.

Thanks

[Mathieu NOE]

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]

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.

Answer 2

[Rajiv Singh]

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]

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.