The function SSICOV.m aims to automatically identify the eigenfrequencies, mode shapes and damping ratios of a line-like structure using ambient vibrations only. The covariance-driven stochastic subspace identification method (SSI-COV) 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. . It has been applied for ambient vibration monitoring of the Lysefjord Bridge  and was compared to the frequency domain decomposition technique . Finally, the algorithm was found accurate enough to visualise the evolution of the bridge eigenfrequencies with the temperature .
The submission file contains:
- A data file BridgeData.mat
- A Matlab Live Script Example1.mlx that illustrates the application of the algorithm.
- A Matlab Live Script Example1_noToolbox.mlx that reproduce Example1 but using the function SSICOV_noToolbox.
- The function SSICOV which is the automated SSI-COV algorithm.
- The function SSICOV_noToolbox which is the automated SSI-COV algorithm but does not use the Statistics and Machine Learning Toolbox. The Linkage algorithm is replaced by the function "PHA_Clustering" by  and the function "cluster" is replaced by the function "Cluster2", which is derived from .
- The function plotStabDiag.m, which plot the stabilization diagram.
Any question, suggestion or comment is welcomed.
 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), 316-329.
 Cheynet, E., Jakobsen, J. B., & Snæbjörnsson, J. (2016).Buffeting response of a suspension bridge in complex terrain. Engineering Structures, 128, 474-487.
 Cheynet, E., Jakobsen, J. B., & Snæbjörnsson, J. (2017).Damping estimation of large wind-sensitive structures.Procedia Engineering, 199, 2047-2053.
 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. 87-93). Springer.
 Yonggang (2021). Fast hierarchical clustering method - PHA (https://www.mathworks.com/matlabcentral/fileexchange/46134-fast-hierarchical-clustering-method-pha), MATLAB Central File Exchange. Retrieved February 4, 2021.
 Eric Ogier (2021). Hierarchical clustering (https://www.mathworks.com/matlabcentral/fileexchange/56844-hierarchical-clustering), MATLAB Central File Exchange. Retrieved February 4, 2021.
Cheynet, E. Operational Modal Analysis with Automated SSI-COV Algorithm. Zenodo, 2020, doi:10.5281/ZENODO.3774061.
very cool stuff! Thanks for sharing.
Thank you very much
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.
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 time-consuming 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 time-steps. If you choose Ts = 10 sec, i.e. 6000 time steps, the function will struggle to compute the Toeplitz matrix while providing not-so-good data for the identification algorithm.
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.
If you have both long records and a sampling frequency of 600 Hz, the cross-covariance 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, 10-20 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.
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 memory-related error. Can you tell me how to solve it?
since I am new in OMA, I'd like to learn this algorithm line by line,
in your modalID function (lines 335-339),
determining the fn, zeta, phi, from eig of A
what is the reference for these?
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 line-like structure must be recorded in at least two different positions.
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));
Requires "linkage" and "cluster" from the Statistics and Machine Learning Toolbox.
The best method I have used ever.
Sometimes measured data is recorded as a column vector. How about putting these lines on line 82?
if Nyy > N
y = y';
thanks for replying .sir i am using matlab 2015 version.
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.
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);
Thank you for the answer
the code is work, if i reduce sampling frequency about 100 times (Fs=1e-4)
By setting a maximal lag of 0.1 sec with a time step of 1e-6 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.
I have a signal with length 0.1 sec and time step 1e-6.
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.
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 non-normalized 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".
That's a good remark! I will try to see if I can implement more rigorously the MAC number for pos = fn0(:)+1-MAC0(:); 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'.
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?
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 "1-MAC0(:)" 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 1-MAC0 is almost zero.
In the paper from Magalhaes et al , the variable pos is indeed without dimension. In the updated version, I have used the definition used in Ref. . 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".
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(:)+1-MAC0(:);
Magalhaes in paper  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
Hi Xinzhe Yuan,
The error is probably triggered because the cross-covariance 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.
Hi there, I got an error:
Undefined function or variable 'T1'.
Error in SSICOV/blockHankel (line 207)
Much appreicated if you can have a look at it.
Find the treasures in MATLAB Central and discover how the community can help you!Start Hunting!