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

Question

Yannick CAIRE 2022-5-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1723590-data-analysis-from-measurements-prbs-injected-to-a-system-how-to-plot-the-bode-of-the-system-res

评论： Yannick CAIRE 2022-10-16

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 个评论
显示无隐藏无

Star Strider 2022-5-20

System Identification Toolbox

Yannick CAIRE 2022-5-20

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2022-5-20

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1723590-data-analysis-from-measurements-prbs-injected-to-a-system-how-to-plot-the-bode-of-the-system-res#answer_968260

在 MATLAB Online 中打开

beam_experiment.mat

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 个评论
显示 1更早的评论隐藏 1更早的评论

Mathieu NOE 2022-5-23

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 2022-5-24

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 2022-5-31

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1723590-data-analysis-from-measurements-prbs-injected-to-a-system-how-to-plot-the-bode-of-the-system-res#answer_974690

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Yannick CAIRE 2022-10-16

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.

请先登录，再进行评论。

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

2 个评论
显示无隐藏无

采纳的回答

3 个评论
显示 1更早的评论隐藏 1更早的评论

更多回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

2 个评论 显示 无隐藏 无

采纳的回答

3 个评论 显示 1更早的评论隐藏 1更早的评论

更多回答（1 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

2 个评论
显示无隐藏无

3 个评论
显示 1更早的评论隐藏 1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论