The function SSICOV.m aims to automatically identify the eigenfrequencies, mode shapes and damping ratios of a linelike structure using ambient vibrations only. The covariancedriven stochastic subspace identification method (SSICOV) is used in combination with a clustering algorithm to automatically analyse the stabilization diagrams.
The algorithm is inspired by the one used by Magalhaes et al. [1]. It has been applied for ambient vibration monitoring of the Lysefjord Bridge [2] and was compared to the frequency domain decomposition technique [3]. Finally, the algorithm was found accurate enough to visualise the evolution of the bridge eigenfrequencies with the temperature [4].
The submission file contains:
 A data file BridgeData.mat
 A Matlab Live Script Example1.mlx that illustrates the application of the algorithm.
 The function SSICOV which is the automated SSICOV algorithm.
 The function plotStabDiag.m, which plot the stabilization diagram.
Any question, suggestion or comment is welcomed.
References
[1] Magalhaes, F., Cunha, A., & Caetano, E. (2009). Online automatic identification of the modal parameters of a long span arch bridge. Mechanical Systems and Signal Processing, 23(2), 316329.
[2] Cheynet, E., Jakobsen, J. B., & Snæbjörnsson, J. (2016).Buffeting response of a suspension bridge in complex terrain. Engineering Structures, 128, 474487.
[3] Cheynet, E., Jakobsen, J. B., & Snæbjörnsson, J. (2017).Damping estimation of large windsensitive structures.Procedia Engineering, 199, 20472053.
[4] Cheynet, E., Snæbjörnsson, J., & Jakobsen, J. B. (2017).Temperature Effects on the Modal Properties of a Suspension Bridge.In Dynamics of Civil Structures, Volume 2 (pp. 8793). Springer.
E. Cheynet (2020). Operational modal analysis with automated SSICOV algorithm (https://github.com/ECheynet/SSICOV/releases/tag/v2.4), GitHub. Retrieved .
Cheynet, E. Operational Modal Analysis with Automated SSICOV Algorithm. Zenodo, 2020, doi:10.5281/ZENODO.3774061.
2.4  See release notes for this release on GitHub: https://github.com/ECheynet/SSICOV/releases/tag/v2.4 

2.3  See release notes for this release on GitHub: https://github.com/ECheynet/SSICOV/releases/tag/v2.3 

2.2  Added Github repository 

2.1  typos corrected + minor modification to improve the performances 

2.0.2  typo 

2.0.1  Set the default value of Ts as 500*dt instead of 20 s to avoid the crashing problem of Matlab if a very high sampling frequency is used. Ideally, Ts should be between two and six times the value of the lowest eigenfrequency of the system. 

2.0.0  The distance parameter "pos" for the clustering algorithm is now properly defined using both the difference in terms of frequencies and the MAC number (thank Mihhail Samusev for the help!). Finally, the description of the submission is updated. 

1.0.4  Updated the definition for the variable "pos" (Cluster analysis) + The default value for the time lag (crosscovariance) is more robustly defined + the variable T1 is preallocated 

1.0.3  Description 

1.0.2  Updated description 

1.0.1  typo 
Create scripts with code, output, and formatted text in a single executable document.
Anna Banas (view profile)
Thank you very much. I didn't notice it at all. Now the program works great. Thank you for the quick answer and for the really great program  it helped me a lot.
E. Cheynet (view profile)
Hi Anna,
It is possible to reduce the maximum lag in the covariance matrix using the optional parameter "Ts". The construction of the Toeplitz matrix on line 203 will be too timeconsuming if the time lag is too large. A large time lag is not suitable for identification purposes either. You can customize the number Ts. A rule of thumb is that Ts can be approximatively equal to 3 x the longest eigen period of the structure. Assuming the first eigenfrequency is 10 Hz, that means an eigen period of 0.1 sec. So Ts = 3x0.1 = 0.3 sec, which is 180 timesteps. If you choose Ts = 10 sec, i.e. 6000 time steps, the function will struggle to compute the Toeplitz matrix while providing notsogood data for the identification algorithm.
Anna Banas (view profile)
Thank you for your response. I thought so, but I did a few tests and unfortunately, this relationship turned out to be not so simple for me. The total number of points of my signal is 15631 and the sampling rate is 600 Hz (26s). For much longer signals and with a lower sampling rate, and I think bigger matrixes the calculations are very fast. Even in the test example, the number of signal samples is much greater than mine. The frequency that is important is around 40 Hz. I can signal downsampled to 120 Hz, but maybe I'm doing something wrong.
E. Cheynet (view profile)
Hi Anna,
If you have both long records and a sampling frequency of 600 Hz, the crosscovariance matrix that is created within the function SSICOV will be too big to be handled in Matlab. If the structure you study has low vibrational frequencies, for example, around 1 Hz, you can reduce considerably the sampling frequency, for example, 1020 Hz and still use long record duration. If your structure has vibrational frequencies around 50 Hz or 100 Hz, it may be good enough to have short record durations, for example, a couple of minutes only.
Anna Banas (view profile)
Hello, I like your program very much and it is very useful. However, I have a problem with one thing. My data is not very long, but when I set my sampling to 600Hz, it takes a long time to calculate and then I get a memoryrelated error. Can you tell me how to solve it?
Benu (view profile)
since I am new in OMA, I'd like to learn this algorithm line by line,
in your modalID function (lines 335339),
determining the fn, zeta, phi, from eig of A
what is the reference for these?
Benu (view profile)
E. Cheynet (view profile)
Note: The problem with Itachi was solved in private discussion. The main conclusions were: (1) the input must be the dynamic acceleration or displacement response, not the free decay response (2) The response along the linelike structure must be recorded in at least two different positions.
Itachi (view profile)
I get this error, I have an 8 second signal with dt= 0.005. Please, how can I get rid of this problem?
Error using SSICOV/blockHankel (line 196)
the IRF must be a 3D matrix with dimensions <M x M x N>
Error in SSICOV (line 95)
[U,S,V] = blockHankel(IRF);
Error in ex (line 59)
[fn0,zeta0,phi0,paraPlot] = SSICOV(Y,t(1));
Mohammad Farshchin (view profile)
Great job!
Modal Analyst (view profile)
Requires "linkage" and "cluster" from the Statistics and Machine Learning Toolbox.
Alex (view profile)
The best method I have used ever.
Kim (view profile)
Sometimes measured data is recorded as a column vector. How about putting these lines on line 82?
if Nyy > N
y = y';
[Nyy,N]= size(y);
end
mati ullah (view profile)
thanks for replying .sir i am using matlab 2015 version.
E. Cheynet (view profile)
Hi mati ullah,
I cannot reproduce the error you get. The compilation is running without crash in Matlab 2018b. I do not know which Matlab version you use but a too old version may be one possible reason for the error message.
mati ullah (view profile)
i have run the example1.mlx with the default data and got the following errors, E. Cheynet sir kindly help me to eliminate these errors as i am new to matlab.
Error in SSICOV (line 149)
[fn3,zeta3,phi3] = myClusterFun(fnS,zetaS,phiS,MACS);
Error in Example1 (line 34)
[fn0,zeta0,phi0,paraPlot] = SSICOV(rz,dt,'Nmin',2,'Nmax',50,'Ts',15);
mmazloumi (view profile)
Thank you for the answer
the code is work, if i reduce sampling frequency about 100 times (Fs=1e4)
E. Cheynet (view profile)
Hi mmazloumi,
By setting a maximal lag of 0.1 sec with a time step of 1e6 sec, you create a covariance matrix that is [ M x M x 1e5], where M is the number of sensors. The resulting Block Toeplitz matrix is way too large to be handled by Matlab. To avoid crashing Matlab, you need to have a number much lower than 1e5.
mmazloumi (view profile)
I have a signal with length 0.1 sec and time step 1e6.
when i put Ts=0.1 sec then there is a error:
"Error using zeros
Out of memory. Type HELP MEMORY for your options.
Error in SSICOV/blockHankel (line 205)
T1 = zeros(M*N1);
Error in SSICOV (line 97)
[U,S,V] = blockHankel(IRF);"
and when i put Ts=1 sec or greater, there is another error:
"Index exceeds matrix dimensions.
Error in SSICOV/NExT (line 270)
IRF(oo,jj,:) = h0(1:M);
Error in SSICOV (line 95)
[IRF,~] = NExT(y,dt,p.Results.Ts,p.Results.methodCOV);"
What is the problem and how i can get true results.
Regards
Wentao Zhao (view profile)
Mihhail Samusev (view profile)
I see, thanks for explanation,
In the article "Cardoso R., Cury A., Barbosa F., 2015  On the advances of automatic modal identification for SHM" they attempt a similar nonnormalized distance between mode estimates, but now including MAC that it scaled with arbitrary constant with [Hz] units, maybe you can take a look at that for inspiration. Maybe, it would be cool for the average user to stick to one of the chosen distance metrics, but for the advanced user to be able to choose among a few or even better, load his own distance function as an optional parameter. Similar to MATLAB's own implementation of "pdist".
E. Cheynet (view profile)
Hei Mihhail,
That's a good remark! I will try to see if I can implement more rigorously the MAC number for pos = fn0(:)+1MAC0(:); in a further version of the submission.
If the MAC number is not included, the clusters are selected based on the criterion 'eps_cluster'. I do not know until which point the cluster analysis will work properly without using the MAC number as an additional criterion. So far, I have been able to distinguish modes separated by only 0.009 Hz for a 1.3 km long span bridge using a low value for 'eps_cluster'.
Mihhail Samusev (view profile)
Thank you for the answer,
Another question, if MAC is not used in the distance, is it then possible to detect closely spaced modes. I see that due to their closeness in frequency they will be clustered together and MAC could be useful to avoid it?
E. Cheynet (view profile)
Hi Mihhail,
That is a very good point! I initially defined the variable "pos" as a quantity with the same dimension as a frequency because it was more intuitive to use the criterion "eps_cluster" this way. I had forgotten to remove the "1MAC0(:)" in the previous version of the function SSICOV. The function still worked well because the dataset investigated in the nested function corresponds already to stables poles, such that 1MAC0 is almost zero.
In the paper from Magalhaes et al [1], the variable pos is indeed without dimension. In the updated version, I have used the definition used in Ref. [1]. However, because the relative distance is now used for the cluster analysis, the criterion "eps_cluster" is much larger in the updated function. I think that the choice of pos = fn0(:) is still acceptable since the MAC number has a limited influence on the cluster analysis. Another argument for using "pos = fn0(:)" is to avoid introducing a dependency between the criteria "eps_cluster" and "eps_MAC".
Mihhail Samusev (view profile)
Hi, thank you very much for your contribution,
I have a question regarding line 450 in SSICOV function,
when you cluster, you firstly calculate positions with:
pos = fn0(:)+1MAC0(:);
Magalhaes in paper [1] uses the same distance, but the frequency is normalized, so is it a mistake in your implementation or you meant to use this distance?
Anyways, really appreciate your work, it helps to develop Structural Health Monitoring toolbox for my university
https://se.mathworks.com/matlabcentral/fileexchange/68988shmtoolbox
E. Cheynet (view profile)
Hi Xinzhe Yuan,
The error is probably triggered because the crosscovariance matrix is not built properly. I will update the function to provide a more robust definition of T1. However, it will probably not solve your problem, but simply trigger a new error. I suggest you to check whether the time lag you use as input (see the optional arguments) is appropriate.
Itachi (view profile)
Xinzhe Yuan (view profile)
Hi there, I got an error:
Undefined function or variable 'T1'.
Error in SSICOV/blockHankel (line 207)
if or(any(isinf(T1(:))),any(isnan(T1(:)))),
Much appreicated if you can have a look at it.
Thiago (view profile)