Error using iddata, Each experiment in the data set must have the same number of input and output measurements (rows)

Hi there,

I am trying to get the vibration frequency with time and acceleration data by system identification using ERA method, when I run the data, there was an error as: error using iddata, Each experiment in the data set must have the same number of input and output measurements (rows). File and code are attached.

Thanks

load C7E_1.mat t acc dt; 

Error using load
Unable to find file or directory 'C7E_1.mat'.

rng('default');

% Add noise
% we add to the excitation and response signals
t = t + randn(size(acc))/1250;
acc = acc + randn(size(acc))/1e11;

% Define outputs from the above function
t = (0:size(t,1)-1)/dt';
X1 = 1e2*t;
Y1 = 1e2*acc;
X0 = X1;
Y0 = Y1(:,1);

subplot(2,1,1)
plot(t,X0)
xlabel('Time (in seconds)')
ylabel('Force (in Newtons)')
grid on;
title('Excitation and Response for the vibration frequency')
subplot(2,1,2)
plot(t,Y0)
xlabel('Time (in seconds)')
ylabel('Displacement (in meters)')
grid on;

%%
% time in a given sample
Ts = 1/t';
% estimation data
% requires System Identification Toolbox.
data_estimation = iddata(Y0(1:1000), X0(1:1000), 1/dt);
% validation data
data_validation = iddata(Y0(1001:2000), X0(1001:2000), 1/dt);


figure
plot(data_estimation)

[~,inputDelay] = max(data_estimation.InputData);
data_estimation = data_estimation(inputDelay:end);

L = length(data_estimation.InputData);
t = seconds(data_estimation.Tstart + (0:L-1)'*Ts);
y1 = data_estimation.OutputData;
tt = timetable(t,y1);


order = 6;
sys = era(tt, order, 'Feedthrough', true);

compare(data_validation,sys);

sys2 = ssest(data_estimation, order, 'Feedthrough', true);

compare(data_validation,sys,sys2);
legend('Validation data','ERA','SSEST');
%%

% Model Evaluation
[~,f] = modalfrf(sys);
[fd, zeta] = modalfit(sys,f,3);

zeta;
fd;


[frf,f] = modalfrf(sys2);
[fd, zeta] = modalfit(sys2,f,3);

%%
%Model Comparison
% lower order
order = 4;
sys = era(tt, order, 'Feedthrough', true);
sys2 = ssest(data_estimation, order, 'Feedthrough', true);

%model comparison with higher order = 6
compare(data_validation,sys,sys2);
legend('Validation data','ERA','SSEST');

[LJsAN]

[Mathieu NOE]

hello

there is something missing in your data : it has only a time and vibration data ; where is the force data ??

then in your code I don't understand what your are doing , taking t for the force data (t is time)

this seems obvious when you do your plot :

[LJsAN]

Hi Mathieu NOE,

Thanks for your reply, I am kinda confused as well because I am new of Matlab only several training but nothing related with ERA method, so I combine everything I know but it seems that not work at all. I am just trying to find the vibration frequency via system identification with time and acceleration data that I had but it seems won't work

[Mathieu NOE]

ok

so what we can at least do here, is to look at the fft spectrum of the vibration data but it tells us the same thing as the time signal itself : it's simply a constant acceleration with a tiny amount of vibration at 6.6 Hz plus higher harmonics

do you believe there was a "real" force excitation or are we simply looking at the background noise here ?

% Fs=2000; % Sampling rate
% N=5000; % Number of data points
% t=[0:1:N-1]/Fs; 
% y=sin(250*2*pi*t)+chirp(t,100,2,500,'quadratic');

load C7E_1.mat; 

Fs = 1/dt;
y = acc;

nfft = 512; 
overlap = 0.75; 
spectrogram_dB_scale = 100; % display this dB range on y axis

[S,F,T] = myspecgram(y, Fs, nfft, overlap); % overlap is 75% of nfft here
sg_dB = 20*log10(abs(S)); % expressed now in dB

% saturation of the dB range :  the lowest displayed level is spectrogram_dB_scale dB below the max level
min_disp_dB = round(max(max(sg_dB))) - spectrogram_dB_scale;
sg_dB(sg_dB<min_disp_dB) = min_disp_dB;

% plot

figure(1); 
plot(t,y) 
xlabel('Time (s)') 
ylabel('Amplitude') 
title('Signal (Acc)')


figure(2); 
subplot(121),
imagesc(T,F,sg_dB) 
hcb = colorbar('vert');
set(get(hcb,'Title'),'String','dB')
set(gca,'YDir','Normal') 
xlabel('Time (secs)') 
ylabel('Freq (Hz)') 
title('Specgram')
colormap('jet');

subplot(122),
sdB = mean(sg_dB,2); % spectrum time averaged = mean along dim 2
plot(sdB,F) 
xlabel('Mean Amplitude (dB)') 
ylabel('Freq (Hz)') 
title('Spectrum')



function  [fft_specgram,freq_vector,time] = myspecgram(signal, Fs, nfft, Overlap)
% FFT peak spectrogram of signal  (example sinus amplitude 1   = 0 dB after fft).
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_specgram = [];
    for ci=1:spectnum
        start = (ci-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_specgram = [fft_specgram abs(fft(sw))*4/nfft];     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_specgram = fft_specgram(select,:);
freq_vector = (select - 1)*Fs/nfft;

% time vector 
% time stamps are defined in the middle of the buffer
time = ((0:spectnum-1)*offset + round(nfft/2))/Fs;

end

[LJsAN]

Thank you very much sir, honestly, I don't think there was a "real" force excitation, I think we are simply looking at the background noise here ? Maybe cuz the acceleration is unfiltered.

are we able to do the identification after we have the ‘original' frequency here?

[LJsAN]

Here is all the data

Thank you

[Mathieu NOE]

hello again

most of the channels will display only a bit of (sensor) noise but a few show indeed interesting results , like C7E_4.mat

I changed a bit my code so that the right plot gives the max amplitude of the spectrogram, so now you see better the "broad" peaks corresponding to the excited eigenmodes

the narrow horizontal lines in the spectrogram plot (on the left) are not relevant, they show some stationnary fixed frequency perturbation (pump or fan or ??) and this kind of narrow band signals are not interresting in modal analysis as they do not cover all frequencies we need to find the eigenmodes of a structure

plot for C7E_4.mat

code for automatic plotting of all mat files :

% use the structure returned by DIR:
P = pwd; % working directory 
S = dir(fullfile(P, 'C*.mat'));

for k = 1:numel(S)

     filename = S(k).name;   % fyi
    F = fullfile(S(k).folder, S(k).name);
    load(F);
    
    % main code 
        Fs = 1/dt;
        y = detrend(acc,'linear');

        %% decimate (if needed)
        % NB : decim = 1 will do nothing (output = input)
        decim = 2;
        if decim>1
            y  = decimate(y ,decim);
            Fs = Fs/decim;
        end

        nfft = 256; 
        overlap = 0.75; 
        spectrogram_dB_scale = 40; % display this dB range on y axis

        [SP,F,T] = myspecgram(y, Fs, nfft, overlap); % overlap is 75% of nfft here
        ind = find(F<=Fs/2.56);
        F = F(ind);
        SP = SP(ind,:);        
        sg_dB = 20*log10(abs(SP)); % expressed now in dB

        % saturation of the dB range :  the lowest displayed level is spectrogram_dB_scale dB below the max level
        min_disp_dB = round(max(max(sg_dB))) - spectrogram_dB_scale;
        sg_dB(sg_dB<min_disp_dB) = min_disp_dB;

        figure(k); 
        subplot(121),
        imagesc(T,F,sg_dB) 
        hcb = colorbar('vert');
        set(get(hcb,'Title'),'String','dB')
        set(gca,'YDir','Normal') 
        xlabel('Time (secs)') 
        ylabel('Freq (Hz)') 
        title(['Specgram for file : ' strrep(filename,'_','-')])
        colormap('jet');

        subplot(122),
%         sdB = mean(sg_dB,2); % spectrum time averaged = mean along dim 2
        sdB = max(sg_dB,[],2); % spectrum peak averaged = max along dim 2
        plot(sdB,F) 
%         xlabel('Mean Amplitude (dB)') 
        xlabel('Max Amplitude (dB)') 
        ylabel('Freq (Hz)') 
        title('Spectrum')

end




%%%%%%%%%%%%%%%%%%%%%%%%
function  [fft_specgram,freq_vector,time] = myspecgram(signal, Fs, nfft, Overlap)
% FFT peak spectrogram of signal  (example sinus amplitude 1   = 0 dB after fft).
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_specgram = [];
    for ci=1:spectnum
        start = (ci-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_specgram = [fft_specgram abs(fft(sw))*4/nfft];     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_specgram = fft_specgram(select,:);
freq_vector = (select - 1)*Fs/nfft;

% time vector 
% time stamps are defined in the middle of the buffer
time = ((0:spectnum-1)*offset + round(nfft/2))/Fs;

end

[LJsAN]

Thank you very much sir, can we just pick some of them like CA7 and C6C to do just an estamite vibration frequency via system identification?

[LJsAN]

The sensors are arranged in two parallel rows on the bridge deck, from left to right, the first row's order is: AFE→C7E→C83→C6F→C80, the second row's order from left to right is: C6C→C6E→CA7→C76→C7F.

C83 and CA7 are located in the middle of the span, AFE and C6C are located at the quarter point.

[Mathieu NOE]

hello again

I don't have the AFE signal, can you provide the data file ?

for several positions (sensors) you have posted multiple files , which ones are the good ones ?

C6C_1.mat

C6C_2.mat

C6C_3.mat

C6C_4.mat

C6C_5.mat

C6C_M1.mat

C6C_M2.mat

C6E_1.mat

C6E_2.mat

C6E_3.mat

C6E_4.mat

C6F_1.mat

C6F_2.mat

C6F_3.mat

C7E_2.mat

C7E_3.mat

C7E_4.mat

C7E_M1.mat

C7E_M2.mat

C80_1.mat

C83_1.mat

CA7_1.mat

CA7_2.mat

CA7_3.mat

[LJsAN]

Hi sir, I don't have AFE signal either, for the sensors data, the following sensors can be considered as good ones

C6C_M1.mat,

C6E_1.mat,

C6F_1.mat,

C7E_M1.mat,

C80_1.mat,

CA7_1.mat

Thank you

[Mathieu NOE]

As far as I can tell / remember , the operationnal modal analysis codes are based on a simple "string" model , in other words we consider only 1D aligned sensors / signals ; I haven't seen any matlab code for a 2D layout

so we are limited to take only one row of sensors at a time

ok so combining your last message infos with previous one :

the first row's order from left to right is : AFE→C7E→C83→C6F→C80,

the second row's order from left to right is: C6C→C6E→CA7→C76→C7F.

We finally can only use fewer points

C6C_M1 , C6E_1 , C6F_1 , C7E_M1 , C80_1 ,CA7_1

first row's : only 3 points valid : C7E_M1→C6F_1→C80_1

second row's: only 3 points valid C6C_M1→C6E_1→CA7_1

[Mathieu NOE]

I did a few checks on the listed "valid" sensors

the problem I have is why the data have different durations ??

we need all data being recorded synchro and same duration otherwise no modal analysis code will work

list = readlines('valid.txt');

for k = 1:numel(list)
    filename = char(list(k));
    F = fullfile(pwd, filename);
    load(F);
    
    % main code 
        Fs = 1/dt;
        acc = detrend(acc,'linear');

        %% decimate (if needed)
        % NB : decim = 1 will do nothing (output = input)
        decim = 2;
        if decim>1
            acc  = decimate(acc ,decim);
            Fs = Fs/decim;
        end
        samples = numel(acc);
        t = (0:samples-1)*dt;
        
        % convert to displacement 
        fc1 = 0.5;
        fc2 = 2;
        [B,A] = butter(2,2*[fc1  fc2]/Fs);
        depl = filter(B,A,acc);
        
        figure(k)
        subplot(211)
        plot(t,depl)
        xlabel('time (s)')
        ylabel(' r_z (m)')
        title(['Displ record from sensor no ' strrep(filename,'_','-')])
        subplot(212)
        pwelch(depl,[],[],[],Fs)
        set(gca,'xscale','log')
        title(['PSD estimate from sensor no ' strrep(filename,'_','-')])
        
end