how can I determine the system's rank so that I can find the transfer function

If I have measured data from a factory (real data), and this data is two signals, the first is an input signal and the second is an output signal.

My question here is how can I determine the system's rank so that I can find the transfer function?

[Markus M.]

What do you mean by rank? The order of the equivalent state-space system or transfer function?

You can try to do system identification and fit a model to your data. You can compute the FFT's of your input and output signals and compute and estimation of your transfer function.

Answer 1

[Mathieu NOE]

hello

if you have a SISO system, there is no "rank" (this applies to MIMO systems) , or rank = 1 by default , unless your output is completely decorreleted from your input (and in this case we are loosing time here)

if you want to plot a Bode diagram from the time data you collected , try this code below

either you have the Signal Processing Toolbox and you can use tfestimate or I offer an alternative as a function in the bottom section

hope it helps

the mat file is attached

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 percent 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); 
% Txy = transfer function (complex), Cxy = coherence, F = freq vector

% Bode plots
figure(1),
subplot(3,1,1),plot(F,20*log10(abs(Txy)));
ylabel('Mag (dB)');
subplot(3,1,2),plot(F,180/pi*(angle(Txy)));
ylabel('Phase (°)');
subplot(3,1,3),plot(F,Cxy);
xlabel('Frequency (Hz)');
ylabel('Coh');

%%%%%%%%%%%%%%%%%%%%%%%
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

Answer 2

[Ahmed Zrigan]

THANK YOU Sir

[Mathieu NOE]

As always, my pleasure !

[Mathieu NOE]

hello again

would you mind accepting my answer ?

tx