A frequency identification method suitable for all-phase spectral refinement of architectural structures

Abstract

Short-time Fourier transformation (STFT) has been widely recognized as an intuitive time-frequency analysis method. However, its application in building structures is constrained by its low accuracy of low-frequency recognition in short data. Accordingly, an all-phase chirp-z transformation (AP-CZT) proposed for frequency recognition in building structures. Specifically, the autoregressive (AR) model, which has a proposed order determination algorithm, is used to solve the problem of data length limitation in the all-phase data process. Subsequently, an algorithm based on absolute inverse proportional function (AIP) is developed for post-refinement frequency correction. To verify its actual application, the ap-CZT method is used to analyze simulating finite element model and white noise feedback data from the actual shake table. The ap-CZT method is proven to be capable of correctly finding high-order frequencies. Moreover, it can accurately identify signal frequencies in short data. Therefore, the ap-CZT method can be applied as a frequency identification method in STFT in the field of building structures.

Keywords:
frequency identification; all-phase; data extrapolation; frequency correction; architectural structures

Graphical Abstract

1 INTRODUCTION

Damaging process of architectural structures has always been an area of great concern for academics and engineers. In general, architectural structures are analyzed through numerical simulation and practical experimentation, with those analyses based on dynamic time-history loading being more accurate. Meanwhile, the analysis of time-history results is still based on the state characteristics of time-domain results, such as the displacement response, acceleration response and force response. These give insufficient description of structural damage process, since they are merely the final presentation of structural state. For frequency domain analysis of structural time-history results, the fast Fourier transform (FFT) is often used. Structural natural frequency and damping constitute the main research content of frequency spectra. Frequency variation before and after experiment, in particular, is analyzed emphatically as the results of shaking table test. Nevertheless, it remains a comparison between various structural states, which is unable to describe the process of structural damage change better (Keliang et al. 2018). Hence, time-frequency analysis (TFA) of architectural structures comes as a vitally effective measure.

TFA has been used widely in the fields like geological survey (Nikoo et al. 2016), seismic exploration (Chenglong 2013; Radad et al. 2015), biomedicine (Sharma et al. 2018), industrial production (Hao et al. 2019) and damage identification (Kunwar et al. 2013). Common basic methods of TFA include (Iatsenko 2015; Sejdić et al. 2018): Short-time Fourier transform (STFT), wavelet variation (WT), Winger-Ville distribution (WV), Hilbert-Huang transform (HHT), etc. Of these, STFT is used extensively in signal processing owing to its simple algorithm, clear principle, easily distinguishable results, and sufficient research. However, its application in the analysis of architectural structures requires overcoming the problems of poor low-frequency resolution.

In reality, STFT is used in the field of architectural structures, and the research related to frequency identification methods remains scarce. The problem is that the first-order natural frequencies are often low for architectural structures, with many buildings showing values below 1 Hz. For the long-period structures (Ding et al. 2014; Fang et al. 2015; Ikeda et al. 2015), the frequencies are only 0.2 Hz (5s) or lower. Scarce researches have been done on the applicability of common frequency identification methods for this bandwidth, especially around 1 Hz. Additionally, achieving good frequency resolution requires a substantial amount of data as support, which is hardly thought achievable in some actual cases. For example, in the shaking table test, the sampling frequency is f _s =1000 Hz. If we want to achieve 0.1 Hz frequency resolution, we need more than 10 seconds’ history data. However, the actual effective time may be less than 5 seconds in the test (Gao et al. 2018). Thus, a high-precision frequency identification analysis method is needed that requires less data.

Initially proposed by Wang and Hou for overcoming the blocking artifacts in image processing, the all-phase digital signal processing mainly addresses the spectrum leakage caused by the truncation error of data, as well as the truncation effects on the truncated signal boundaries (Amiri 1983; Wang et al. 2001). Meanwhile, STFT needs to truncate the original data multiple times. The all-phase FFT (apFFT) spectral analysis is one effective processing method for suppressing such error. To enhance the accuracy of frequency identification, an integrated approach of spectral refinement and frequency correction can be used without increasing the data length. Common methods of spectral refinement include the Zoom-FFT, chirp-z transform(Rabiner et al. 1969) (CZT), FFT-FS, etc (Mao et al. 2012). All of these methods approximately interpolate the frequency values between spectral lines via a certain mathematical approach, thereby improving the pseudo frequency resolution of signals. Excessively high degree of refinement, however, is detrimental to the frequency identification since it results in excessively smooth spectrum. Accordingly, spectral correction methods can be used, to improve the identification accuracy further. A category of direct methods is to correct the frequency using amplitude spectrum, such as the interpolated FFT amplitude spectrum(Jain et al. 1979) and the energy centrobaric correction method (Ding and Jiang 2001). A principal factor affecting their estimation accuracy is the inter-spectrum interference resulting from spectrum leakage (Xiao and Dong 2013). The correction methods based on apFFT have better accuracy owing to their excellent performance in suppressing spectrum leakage.

As the above studies show, the methods of all-phase processing, spectral refinement and frequency correction are interrelated essentially. Among them, all-phase data processing technology is the basis, which suppresses spectrum leakage, improves spectral refinement frequency resolution, and enhances frequency correction accuracy. Meanwhile, spectral refinement limits frequency correction method error; and frequency correction method improves frequency identification accuracy of spectral refinement. Therefore, this paper proposes a frequency identification method (ap-CZT) that incorporates the above three techniques and applies it to time-history analysis of architectural structures. Moreover, some key issues related to methodological construction with the ap-CZT are resolved.

2 Principle and algorithm

In terms of methodological implementation, there are two major problems under resolution: (1). Prior to performing the all-phase processing, the preprocessing of signals for architectural structures is necessary. (2). the refined spectrum should be subjected to the peak frequency correction. To this end, the construction principle of the all-phase preprocessing-based CZT spectral refinement method is described below.

2.1 Discrete sequence preprocessing

The input data resulting from all-phase preprocessing must be required to have a length of 2N-1, while the output data is N in length. In other words, although the subsequent analyses are all transformations based on N number of data, the data size used is 2N-1, implying that the N-1 number of data is not used in the subsequent transformations. For architectural feedback signals, their length is not excessively long, so halving data size is unacceptable. In TFA, the number of data samples determines frequency accuracy, so the entire samples should be used as much as possible. Therefore, a data processing method is needed extending the N number of data to 2N-1. Common data extension approaches include original data reuse, zero padding, and data extrapolation.

From a structural perspective, this method easily judges the situation of extended data after all-phase processing. A sequence of data with 2N-1 length is defined as shown in equation (1).

$x^{2N-1}(i),i=1-N,\dots ,N-1$ (1)

After undergoing the all-phase preprocessing, the data with a sample number of N is $x_{ap}^N\left(m\right),m=\mathrm{0,1,}\dots ,N-1$. Later, in the case of rectangular window preprocessing, the processed data $x_{ap}^N\left(m\right)$ are expressed as

$x_{ap}^N(m)=(N-m)x^{2N-1}(m)+m{x}^{2N-1}(m-N),m=\mathrm{0,1,}\dots ,N-1$ (2)

and, it will become the basic expression in this section.

1. Data shift by abandoning the first bit (DSAF)

The simplest method of data extension is the repeated use of known data, where original data is extended through shift/flipping measures. This approach is equivalent to using the first N data $x^N\left(i\right),i=1-N,\dots \mathrm{,0}$ of $x^{2N-1}\left(i\right)$ as the original data, and abandoning the first bit from the last N-1 data. Its structure is:

$x_{DSAF}^{2N-1}(i)=\{\begin{array}{l}{x}^N(i),i=1-N,\dots \mathrm{,0}\hfill \\ x^N(i-N+1),i=\mathrm{1,}\dots ,N-1\hfill \end{array}$ (3)

Then, the data after all-phase (equation (2)) processing are:

$x_{ap}^N(m)=\{\begin{array}{l}N{x}_{DSAF}^{2N-1}(0),m=0\hfill \\ (N-m)x_{DSAF}^{2N-1}(m-N+1)+m{x}_{DSAF}^{2N-1}(m-N),m=\mathrm{1,}\dots ,N-1\hfill \end{array}$ (4)

Compared to equation (2), the equation (4) basically retains the data structure of all-phase processing as much as possible. However, in essence, it is still based on the principle of superimposing original data in different ways, which is less effective than the all-phase processing of 2N-1 data. Nevertheless, it (DSAF) can be used as a reference method for comparison owing to its simple construction.

2. Zero padding (ZP)

Zero padding is used commonly in the FFT to pad each row of X with zeros so that the length of each row is the next higher power of 2 from the current length. Extending data to a length of 2N-1 only requires adding N-1 zeros before or after the sequence $x^N\left(i\right),i=1-N,\dots \mathrm{,0}$. Later zero padding is considered herein (prior zero padding is similar in results):

$x_{ZP}^{2N-1}(i)=\{\begin{array}{l}{x}^{2N-1}(i),i=1-N,\dots \mathrm{,0}\hfill \\ \mathrm{0,}i=\mathrm{1,}\dots ,N-1\hfill \end{array}$ (5)

Then, the data after all-phase processing are:

$x_{ap}^N(m)=\{\begin{array}{l}N{x}_{ZP}^{2N-1}(0),m=0\hfill \\ m{x}_{ZP}^{2N-1}(m-N),m=\mathrm{1,}\dots ,N-1\hfill \end{array}$ (6)

From the equation (6), it is clear that after all-phase processing, the later zero padding approach is actually equivalent to adding $x_{}^{2N-1}\left(i\right)$ with an eccentric triangular window, and cyclically shifting the last data to the first bit. Zero padding is used as one of the comparison methods, since it is a common data extension approach.

3. AR model-based data extrapolation (AREX)

Data extrapolation allows increase of data length while keeping the sampling frequency unchanged. Accordingly, the data extension problem can be described as follows: given a time-domain sequence $x^N\left(n\right),n=\mathrm{1,2}\dots N$, with an original length of N, its spectrum after Fourier transform is $X^N\left(k\right),k=\mathrm{1,2}\dots N$, and its amplitude spectrum $A^N\left(k\right),k=\mathrm{1,2}\dots N$. Estimate its time-domain sequence ${\widehat{x}}^{2N-1}\left(n\right),n=\mathrm{1,2}\dots 2N-1$ after data extension by N-1.

AR models are models in which the new data are derived based on a linear combination of the previous p data. Accordingly, for a known time-domain sequence$x^N\left(n\right),n=\mathrm{1,2}\dots N$, we can get the estimated sequence ${\widehat{x}}^{2N-1}\left(n\right),n=\mathrm{1,2}\dots 2N-1$ through N-1 times of predictions without self-correction (Zhou et al. 1999). However, since it is a prediction without self-correction, predicting the estimate from the estimated value is unavoidable. Hence, the AR models are essentially smooth, and the frequency domain characteristics of original sequence will be continued into the prediction part. It is thus necessary to ensure that the spectral characteristics of sequence $x^{2N-1}\left(n\right),n=\mathrm{1,2}\dots 2N-1$ are similar throughout the first and the second half parts. This can be guaranteed in the case of short data in this study.

The estimate $\widehat{x}\left(n\right)$ of AR-based prediction models can be expressed as a weighted sum of past p values:

$\widehat{x}(n)=-{\displaystyle \sum _{k=1}^p{a}^p}(k)x(n-k)$ (7)

Where p denotes the order of AR (p) models, and $a^p\left(k\right),k=\mathrm{1,}\dots ,p$ are the model parameters. Among many ways of estimating model parameters, the maximum entropy spectral recursion proposed by J. P. Burg (Burg 1975) is adopted as the basic concept in this study, owing to its advantages like simple calculation, high resolution and adaptability to short time sequences. This paper uses its modified algorithm for the AR model parameter estimation(Shen et al. 2008). Accordingly, the time-domain sequence extension method herein can be described as follows:

• Determine the optimum order p of AR (p) models.

• Using the algorithm in reference(Shen et al. 2008), estimate the AR model parameters$a^p\left(k\right),k=\mathrm{1,}\dots ,p$ by the original time-domain sequence $x^N\left(n\right),n=\mathrm{1,2}\dots N$.

• Initialize the data ${\widehat{x}}^{2N-1}\left(n\right)=x^N\left(n\right),n=\mathrm{1,2}\dots N$.

• Calculate the rest of ${\widehat{x}}^{2N-1}\left(n\right),n=N+\mathrm{1,2}\dots 2N-1$ using equation (8).

${\widehat{x}}^{2N-1}(n)=-{\displaystyle \sum _{k=1}^p{a}^p}(k){\widehat{x}}^{2N-1}(n-k),n=N+\mathrm{1,}\dots \mathrm{,2}N-1$ (8)

In the computational course, the AR model parameter $a^p\left(k\right)$ is constant. Actually, it can be re-estimated along with the predicted data. However, updating $a^p\left(k\right)$ without self-correction may bring large errors. Similarly, in the data extrapolation process, it is also unavoidable to use unverified estimates for data prediction. Therefore, the selection of an appropriate order p is crucial.

Common order determination criteria for AR models, such as Akaike’s Information Criterion (AIC), Bayesian Information Criterion (BIC) and Final Prediction Error (FPE), have varying degrees of adaptability to different problems. The problem herein is to obtain ${\widehat{x}}^{2N-1}\left(n\right)$ through the time-domain sequence $x^N\left(n\right)$, and to make it as similar as possible to ${\widehat{x}}^{2N-1}\left(n\right)$. A similar purpose is to make the spectral characteristics of estimated sequence ${\widehat{x}}^{2N-1}\left(n\right)$ after all-phase processing change as little as possible. Hence, an order determination method suitable for the present problem can be developed from an energy perspective.

In order to make the spectrum as similar as possible, energy should be concentrated on the effective frequency bandwide. Therefore, an energy concentration criterion (ECC) based on FFT amplitude spectrum is constructed, and its expression is shown in equation (9).

$ECC(p)=\frac{Max\left[{\left({\widehat{A}}_p^{2N-1}(m)\right)}^2\right]}{\frac{1}{N}{\displaystyle \sum _{n=k_{\lim }+1}^N{\left[{\widehat{A}}_p^{2N-1}(n)\right]}^2}},m=\mathrm{1,..,}{k}_{\lim }$ (9)

Where, Max[x] represents the maximum value of x. And, $k_{\lim }$ is the limit of effective frequency bandwide. For architectural structures, $f_s=1000Hz$ is a common sampling frequency, and its quarter (125Hz) is enough to cover the main frequency of the structures. So, $k_{\lim }$ can be N/4. Equation (9) can be described as the ratio of the maximum value of the square of the frequency amplitude in the effective frequency range (approximately understood as energy) to the average value of the square of the amplitude of the non-effective frequency range with the change of the order p. What’s more, when the maximum value of ECC is taken, the corresponding p value the optimum order of AR model, since the greater the value, the more energy is concentrated in the effective frequency range.

The following is a comparative analysis of $ECC\left(p\right)$,$FPE\left(p\right)$,$AIC\left(p\right)$, and$BIC\left(p\right)$, and their expressions are defined as follows:

$FPE(p)=\frac{N+p}{N-p}{\widehat{\sigma }}_p^2$ (10)

$AIC(p)=N\ln {\widehat{\sigma }}_p^2+2p$ (11)

$BIC(p)=N\ln {\widehat{\sigma }}_p^2+p\ln N$ (12)

where,${\widehat{\sigma }}_p^2$ is estimation of white noise variance for p-order AR model.

Now, a third-order dominant frequency signal is expressed by equation (13):

$x^{2N-1}(n)=3\sin \left(12.55*2\pi n/f_s\right)+\sin \left(17.55*2\pi n/f_s\right)+2\cos \left(21.75*2\pi n/f_s\right)$ (13)

where $f_s=1000Hz$,$n=\mathrm{1,}\dots \mathrm{,2}N-1$. The first N data $x^N\left(n\right),n=\mathrm{1,}\dots ,N$ of $x^{2N-1}\left(n\right)$ are the original data. By selecting different orders p, that estimation sequence ${\widehat{x}}^{2N-1}\left(n\right)$ in different orders p is obtained according to the AREX algorithm. And, the correlation between the estimated sequence ${\widehat{x}}^{2N-1}\left(n\right)$ and the target sequence $x^{2N-1}\left(n\right)$ is expressed by Pearson correlation coefficient.

$err(p)={\displaystyle \sum _{n=0}^{2N-1}{\left[x^{2N-1}(n)-{\widehat{x}}^{2N-1}(n)\right]}^2}$ (14)

When $err$(p) takes the minimum value, p is the true optimum order. Morever, combined with engineering experience, the range of orders for p is $\left[N/\mathrm{3,}N-1\right]$.

Assuming data length N=128, the order p of ECC method and other three methods are listed in Figure 1, for convenience of comparison, all of them are normalized.

Figure 1
Contrast diagram of optimal order p determined by different determination methods.

It can be observed from the Figure 1 that at a small data length N=128, the values within the order value range are basically unchanged with the three methods of FPE, BIC and AIC, while the actual error (err) values change within the range. At a minimum err value; the optimum order p is 74. Contrastively, with the ECC method in this paper, the order is set to be 72. Despite not being optimal, it is far more accurate than the other three methods.

In reality, signals are always accompanied by noise, so it is necessary to consider the performance of the ECC method under noisy signals. The sequence in equation (13) is subjected to 1000 cycles of Monte Carlo analysis separately using three different signal-to-noise ratio(SNR) additive white noises (SNR = 10 db, SNR = 0 db, SNR = -10 db) by four different order determination methods. Table 1 presents the relevant order determination results.

Table 1
The optimal order p determined by various determination methods.

SNR	ECC	FPE	AIC	BIC	err
10db	59	124	127	127	84
0db	68	126	127	127	66
-10db	78	118	127	127	51

According to Table 1, orders determined by FPE, AIC and BIC methods are all somewhat high, which are unable to truly reflect the optimum model order. In contrast, the order determined by the proposed ECC method is closer to the best order. The error corresponding to the ECC method is incredibly close to the minimum error, while other methods show larger deviations. Thus, the proposed ECC method is applicable to the AR model order determination of extended data in architectural structure signals.

4. Verification of data extrapolation methods

Based on the foregoing analyses, there are three methods that can be used for order extension, namely the data shift by abandoning the first bit (DSAF), the zero-padding (ZP) and the AR model-based data extrapolation (AREX). Using equation (13), the three kinds of data extension methods are examined for verification of respective performances. Supposing that N=128 and N=256, and the AR model order determined via the present ECC method is p=72 and p=258, respectively. Then, the sequence extended by the three methods is ${\widehat{x}}^{2N-1}\left(n\right)$(DSAF, ZP, AREX), and the original sequence is $x^{2N-1}\left(n\right)$(ORI, this means that the second half of the data is obtained by equation (13)). Figure 2 presents their results after rectangular window all-phase data processing.

Figure 2
The waveform results after rectangular window all-phase data processing.

The waveforms of AREX and ORI sequences have some deviations, but they are far less than those of other methods at a data length of N= 128, as shown in Figure 2a. However, at a data length of N= 512, the waveforms of AREX and ORI almost coincide, as shown in Figure 2b.

The estimated sequence ${\widehat{x}}^{2N-1}\left(n\right)$ obtained by different methods and the original sequence $x^{2N-1}\left(n\right)$ are analyzed by apFFT, and the amplitude spectrum is shown in Figure 3

Figure 3
The apFFT amplitude spectrum by three estimated sequences and original sequence.

According to the amplitude spectra, the principal frequencies identified by DSAF and ZP are all shifted at a data length of N= 128, as shown in Figure 3a, while that identified by AREX does not deviate, showing an amplitude almost similar to ORI. At a long data length of N= 512, as shown in Figure 3b, none of the three methods exhibits frequency shift. Nevertheless, as can be seen from the black box part, apFFT suppresses the spectrum leakage with both AREX and ORI, while the remaining methods fail to reflect the advantages of the all-phase processing well. In addition, the amplitude spectra of AREX and ORI almost overlap. Thus, clearly, AREX is a preferable method of data extension prior to all-phase processing.

2.2 Frequency correction method based on axisymmetric function

Correction utilizing relationships between spectral lines is a common measure of spectrum correction. Common such methods, such as the energy centrobaric (Ding and Jiang 2001) and interpolated methods, are based mainly on the DFT spectrum, but the refinement of spectrum has been studied insufficiently. In this section, a frequency correction method based on three-point coordinates can be mainly introduced.

To unify the expression for various methods, the spectral symbols are unified.

Figure 4
The apFFT amplitude spectrum by three estimated sequences and original sequence.

In Figure 4, the coordinates of the main spectrum, the previous and back post points are$K\left[x\left(K\right),y\left(K\right)\right]$,$M\left[x\left(M\right),y\left(M\right)\right]$, and $N\left[x\left(N\right),y\left(N\right)\right]$, respectively. Also, it must satisfy the inequality $y\left(K\right)\ge y\left(M\right)\ge y\left(N\right)$. From the above, the correction frequency can be obtained. Besides, point O is the correction point after frequency correction, and its abscissa $x\left(O\right)$ is the correction frequency.

Defining the frequency resolution of the DFT spectrum$\Delta f$, coordinates of M and N points can be expressed as

$\{\begin{array}{l}y(M)=my(K)\hfill \\ y(N)=ny(K)\hfill \\ x(M)=x(K)\pm \Delta f\hfill \\ x(N)=x(K)\pm \Delta f\hfill \end{array}$ (15)

where $m\in \left[\mathrm{0,1}\right]$, $n\in \left[\mathrm{0,}m\right]$. Because of the inequality$y\left(K\right)\ge y\left(M\right)\ge y\left(N\right)$ , the distance from point M to point K must be greater than or equal to the distance from point N to point K, then the correction frequency $x\left(O\right)$ should be between $x\left(K\right)$ and$x\left(M\right)$.

In order to express more concisely, according to the above definition,$x\left(O\right)$ can be defined as:

$x(O)=x(K)+k\Delta f$ (16)

where, k was defined as a correction scale factor, which can represent the extent to which the correction frequency$x\left(O\right)$ deviates from $x\left(K\right)$.

1. Principle of the frequency correction method

In this section, the correction method can assume that the attenuation spectral line function$f\left(x\right)$ is an axisymmetric function with no inflection point in a frequency resolution pre and post the main spectral line. Then, on the DFT spectrum of the signal, the main frequency spectrum point, the previous secondary frequency spectrum point and the post-secondary frequency spectrum point should be on the axisymmetric function. Thus, by specifying function form and considering known three-point coordinates, specific function expression can be determined. Then, its abscissa of symmetric axis was the corrected frequency value.

Besides, with equation (15),$f\left(x\right)$ is an axisymmetric function without inflection point, so it can be shown as $x\left(O\right)\in \left[x\left(K\right)-\Delta f/\mathrm{2,}x\left(K\right)+\Delta f/2\right]$. Therefore, the$f\left(x\right)$ boundary conditions must be shown as:

$\{\begin{array}{l}x(O)=x(K)\pm 0.5\Delta f,m=1\hfill \\ x(O)=x(K),m=n\hfill \end{array}$ (17)

Absolute value function is a kind of common axisymmetric function. In this paper, we constructe a kind of absolute value function:

$y=a|x-b{|}^{\alpha }+c$ (18)

where, a, b, and c can be calculated by three-point coordinates, and $\alpha$ is the shape factor. For example, when $\alpha =2$, it is quadratic function; when $\alpha =1$, it is primary absolute value function; when $\alpha =-1$, it is absolute value of inverse proportional function(AIP). Combining equation (15)-(18), the expression for the correction scale factor k is:

$k=\mathrm{sgn}(x(M)-x(N))\frac{n-m}{2(n-1)},\alpha =1$ (19)

$k=\mathrm{sgn}(x(M)-x(N))\frac{n-m}{2(m+n-2)},\alpha =2$ (20)

$k=\mathrm{sgn}(x(M)-x(N))\frac{m-n-{[(8+m-9n)(m-n)]}^{1/2}}{4(n-1)},\alpha =-1$ (21)

2. Performance of axisymmetric function-based frequency correction method

First of all, applying equation (15) and (16) to interpolation method and energy centrobaric correction method, k can be shown as following.

For interpolation method with rectangular window that uses three DFT magnitude samples(Jacobsen and Kootsookos 2007), k can be same as equation (20):

$k=\mathrm{sgn}[x(M)-x(N)]\frac{n-m}{2(m+n-2)}$ (22)

Besides, according to energy centrobaric correction methord (Ding and Jiang 2001) by using three DFT magnitude samples, k can be shown:

$k=\mathrm{sgn}[x(M)-x(N)]\frac{m-n}{m+n+1}$ (23)

To analyze the principle of different methods, it needs to understand the trend of k changing with m and n. To analyze the principle of different methods, it needs to understand the trend of k changing with m and n. There are three methods by three-point, which are the interpolation method, energy centrobaric correction method, and the AIP method proposed.

Figure 5
Relationships among k, m and n.

The foremost difference of Figure 5b from Figures 5a and Figures 5c is the boundary value assignment when m=1. This is attributed primarily to the different underlying concepts of the methods, but pros and cons comparison is not the focus of this paper. In addition, it is shown in Figures 5c that the surface protrudes obviously near the m = n line. Thus, with respect to AIP method, when the difference between m and n is small, it will have a big impact on k. However, the change in k value tends to be smooth with the increasing difference between m and n. After CZT refinement, the amplitude curve of spectrum becomes smoother, so that the m and n values are close, which approach 1. In the case of the energy centrobaric correction, as shown in Figure 5b, the k value changes slightly. The interpolation method shown in Figures 5a has some improvements, as it conforms to this paper’s axisymmetric function concept. Further, the AIP correction method has high sensitivity near m = n, which is quite suitable for frequency correction following the CZT spectral refinement.

2.3 All-phase CZT spectral refinement (ap-CZT)

Based on the above analysis, the all-phase data preprocessing procedure can suppress spectrum leakage; the frequency resolution can be improved by the spectrum refinement of CZT, whose accuracy of frequency identification may be further improved by a frequency correction method of AIP. Therefore, we suggest the Ap-CZT method that can be officially described as:

1. A data length extension from N to 2N-1 with AREX can be necessary before all-phase data preprocessing procedure.

2. Carry out all-phase processing on the data(Xiao and Dong 2013), and rectangular window, single-window or double-window can be selected as required.

3. Apply the CZT(Rabiner et al. 1969) spectral refinement algorithm with determining the refined frequency bandwidth $\left[f_1,f_2\right]$ and the refining multiple D.

4. Correct the frequency using the AIP method or other frequency correction method based on axisymmetric function.

5. Output the main frequency results.

3 Results and discussion of experiments

For verification purpose, the proposed ap-CZT algorithm is numerically simulated and applied to actual engineering analysis. Initially, the performance of ap-CZT is compared with apFFT and CZT, in order to reflect its advantages over single algorithms. Next, it applies to handling FEA data. In comparison to the finite element modal method, the present method exhibits a frequency identification error of less than 5%, and allows identification of the hardly identifiable high-order mode frequencies. Finally, the ap-CZT method is applied to analyzing shaking table test data. Under high-noise signals, the conventional FFT is unable to identify the specific frequencies even with large data size. In contrast, the proposed method can achieve accurate frequency identification with small data size.

3.1 Numerical spectrum simulation results with ap-CZT

For architectural structure data, they need to contain at least the third-order mode, and the amplitude and initial phase of third-order principal frequencies should be different. Moreover, the third-order principal frequencies are similar and low, which generally do not exceed 50 Hz. Regarding the acquisition equipment, its sampling frequency is usually between 100~1000 Hz. Accordingly, after considering the characteristics of architectural structure data, a real sequence of third-order principal frequencies is created.

$x(n)={\displaystyle \sum _{i=1}^3{A}_i}\sin \left(f_i\times 2\pi n/f_s+{\phi }_i\right),n=\mathrm{1,}\dots ,N$ (24)

The length of data is N = 1024, and the sampling frequency is $f_s$=1000Hz. The remaining parameters are detailed in equation (25). Noiseless and noisy scenarios are considered, respectively.

$\{\begin{array}{l}{A}_1=4\hfill \\ f_1=0.432\hfill \\ {\phi }_1=0.5\pi \hfill \end{array},\{\begin{array}{l}{A}_2=2\hfill \\ f_2=1.692\hfill \\ {\phi }_2=0.25\pi \hfill \end{array},\{\begin{array}{l}{A}_3=1\hfill \\ f_3=2.445\hfill \\ {\phi }_3=0\hfill \end{array}$ (25)

1. Noiseless data analysis

From the perspective of verification, in order to ensure the consistency of initial input, the following data are adopted. The same data of apFFT and ap-CZT are extrapolated the N number of data to 2N-1 with the AREX method. Besides, the refined frequency bandwidth is$\left[0Hz\mathrm{,8}Hz\right]$ and the refining multiple D is 8. So, their amplitude spectrums can be shown in Figure 6.

Figure 6
The spectrum analysis comparison charts among apFFT, CZT, and ap-CZT, N=1024.

The signal has third-order frequencies of 0.432 Hz, 1.692 Hz and 2.445 Hz. As is clear from figure 6a, the apFFT method is unable to identify the third-order frequency information at a data length of N = 1024. The main constraint of identification accuracy is the insufficient frequency resolution attributed to the data size.

According to Figure 6b, the CZT method can distinguish the second-order frequencies, with errors of 44.68% and 10.82%. However, due to the oscillation attenuation of principal frequencies caused by spectrum leakage, the identification of third-order main frequencies is impacted. Further, the data length will be set to N = 2048 while keeping the other parameters unchanged. In this setting, CZT method can identify all three order frequencies, but false peaks are present between the second and third orders, which will affect the frequency identification severely.

It can be seen from Figure 6c that the proposed ap-CZT method can distinguish the third-order frequencies accurately, without any missing or false frequency. The relevant third-order frequency errors are 4.861%, 1.478% and 3.272%. Respecting the accuracy of frequency quantity identification, the ap-CZT allows accurate identification of quantity at small data size, which is more accurate than the apFFT and CTZ refinement methods. To sum up, ap-CZT integrates the advantages of the two methods remarkably. Respecting the frequency identification accuracy, the third-order frequency errors are all below 5% with the present method, which basically meets the requirements for use in architectural engineering. In first-order frequency identification, its accuracy is much higher than the other two methods. This is mainly owing to the principal frequency correction method of this paper. Comparison with Figure 6 reveals that for the images after 3 Hz, ap-CZT exhibits highly strong frequency suppression in the non-principal frequency bands, which is basically near the zero line. In this way, the energy is concentrated on the principal spectral lines, making the frequency correction more accurate.

2. Noisy data analysis

In this section, the proposed ap-CZT method is verified by adding noise to the signal in equation (24). The signal SNR values are set at -20 db, -10 db, 0 db, 10 db and 20 db, and the identified frequencies of 1000 random noise signals are averaged. Table.3 presents the details.

Table 2
The dominant frequency identified by ap-CZT under different SNR

SNR(db)	-20	-10	0	10	20
1st(0.432Hz)	0.5325	0.5001	0.4980	0.4970	0.4703
2nd(1.692Hz)	1.7471	1.7193	1.6864	1.6884	1.6895
3rd(2.445Hz)	2.8430	2.6251	2.5321	2.5295	2.5284

As shown in Table 2, at positive SNRs, AP-CZT results are within acceptable range. The identified second and third-order frequencies differ little from those in the noiseless scenario. The errors with first-order frequencies are due to the excessively short data length, excessively low first-order frequencies, incomplete waveforms and severe impact of noise. Under negative SNRs, there is no misidentification of main frequency or divergence, despite slightly worse identification accuracy. It is remarkable, though, that due to the limitations of ECC rules in this paper, at SNRs <-20, the results will be largely erroneous, where longer data are needed to lower the impact of error.

3.2 FEA simulation results of architectural structures with ap-CZT

In this section, the performance of the proposed ap-CZT method in practical structure engineering is verified by utilizing the numerical modeling analysis results of a certain industrial plant structure. The plant has a complex structure, with uneven stiffness and mass distributions, so high-order modes are involved in the structural vibrations. Figure 7 illustrates its shaking table model.

Figure 7
The model of industrial plant structure.

As can be seen, it is a scaled model of the shaking table test. Relevant modeling and structural modal computation are accomplished by SAP2000, where the Y-direction modals are considered only. In Figure 8, the highly contributing modals and corresponding frequencies are presented.

Figure 8
The shape and frequency of Y-direction modals.

The initial small deformation is applied to the Y direction of the structure, which is suddenly released, and the structure is free to attenuate vibration. The parameters of apFFT can be defined as follows. Sampling frequency $f_s=1000Hz$, the length of data$N=1024$, the optimal order p=539 with ECC method, the refined frequency bandwidth is$\left[0Hz\mathrm{,8}Hz\right]$ and the refining multiple D is 16. Further, it needs to be recorded that the acceleration time-history data of sampling point in the Y-direction, which can be found in Figure 7a, and its results after apFFT can be shown in Figure 9 (Normalized amplitude spectrum).

Figure 9
The amplitude spectrum of acceleration data in Y-direction.

It can be seen from Figure 9 that the ap-CZT spectral analysis method can accurately distinguish the first four orders of frequencies that contribute to the Y-direction. The errors of these frequencies are 2.45%, 0.73%, 0.28% and 1.22%, respectively, which are in line with the accuracy requirements in the architectural structure field. Therefore, it has a high computational efficiency and can be used for analyzing structural FEA time history data in STFT.

3.3 Shaking table test results with ap-CZT

To verify the performance of the proposed ap-CZT method in actual architectural structures, a frequency analysis is performed on the initial white noise acceleration signal of a three-story steel brace-concrete frame shaking table test model (F-BRC) described in reference [14]. The relevant FEA modal analysis results, which represent the first 4 frequencies of the structure, are 5.52Hz, 21.25Hz, 33.31Hz, and 41.58Hz, and the experimental structural model has been shown in Figure 10. The problem is that its data has high noise intensity, and noises are colored, polluting the identifying frequencies, leading to poor frequency identification, as shown in Figure 11a.

Figure 10
The experimental structural model.

In the shaking table experiment, unidirectional white noise signals are applied to the model structure, and acceleration sensor signals are collected at a data sampling frequency of. With the conventional FFT, all data (N=65536) are used, and Figure 11a displays relevant spectral analysis results (normalized). With the ap-CZT method, only shorter data (N=256) are used. Meanwhile, the order of AR models is p = 103, the frequency band is, and the refinement multiple is D = 16. The relevant results (normalization) are presented in Figure 11b.

Figure 11
The spectrum of white noise acceleration signal in shaking table test.

According to Figure 11a, there are severe signal noise interferences. Except for the barely observable peaks in the first- and fourth-order cycles, the positions of second- and third-order peaks are unclear. Nonetheless, it is still observed that the fourth-order modes prevail. From Figure 11b, it can be seen that the proposed method can distinguish between four mode frequencies, and the peaks of four frequency types are apparent. It is difficult to rule out the construction errors, since the main structure is made of concrete. Thus, the actual frequencies of the structure may somewhat differ from the numerical model calculations, and the actual error comparison is of little importance. Nevertheless, their approximate range is still within the scope of numerical modal frequencies. Furthermore, with an extremely small data use (N=256), ap-CZT is suited to shaking table test analysis of scaled models.

4 CONCLUSIONS

• This paper proposes a spectral analysis method applicable to architectural structure analysis by integrating all-phase data preprocessing with the CZT.

• Taking into consideration the characteristics of all-phase data processing, this paper compares the common data extension approaches, and concludes that the AR model-based data extension (AREX) best reflects the advantages of all-phase processing technology. Meanwhile, a novel order-defining method (ECC) suitable for data extension is being developed.

• To address the inefficient correction with the energy centrobaric and Interpolation methods after spectral refinement, a frequency correction method (AIP) following CZT refinement is proposed based on the idea of axisymmetric function.

• Verification with both numerical and actual data shows that the proposed ap-CZT is a frequency identification method for architectural structures. Featuring high accuracy at short data lengths, it can be used as a frequency identification method of STFT, which requires further study.

References

• Amiri WZ-HH (1983) THE RECONSTRUCTION OF SUB-NYQUIST SAMPLING PAL SIGNAL BY TWO DIMENSIONAL ALIASING DIGITAL FILTER. J Tianjin Univ 1:
• Burg J (1975) Maximum Entropy Spectral Analysis. Ph Thesis Stanf Univ
• Chenglong PRLBS (2013) Review on time-frequency analysis technique and its application in seismic exploration. Lithol Reserv 25:92-96
• Ding JM, Wu HL, Zhao X (2014) Current situation and discussion of structural design for super high-rise buildings above 250 m in China. J Build Struct 35:1-7
• Ding K, Jiang L (2001) Energy centrobaric correction method for discrete spectrum. J Vib Eng 14:354-358
• Fang MX, Yang ZY, Zhi LH (2015) Inverse analysis of across-wind loads on super tall buildings. J Vib Shock 34:35-41
• Gao X, Jianqin LI, Liu C, Yanglong LI (2018) Mechanism of dynamic inelastic torsion and anti-tosional design strategy of steel-braced concrete frames. J Build Struct 39:44-53. https://doi.org/10.14006/j.jzjgxb.2018.02.005
  » https://doi.org/10.14006/j.jzjgxb.2018.02.005
• Hao G, Li F, Bai Y, Wang W (2019) Time-frequency Analysis of Non-Stationary Signal Based on NDSST. Geomat. Inf. Sci. Wuhan Univ. 44:941-948
• Iatsenko D (2015) Linear Time-Frequency Analysis. In: Iatsenko D (ed) Nonlinear Mode Decomposition: Theory and Applications. Springer International Publishing, Cham, pp 7-42
• Ikeda A, Fujita K, Takewaki I (2015) Reliability of System Identification Technique in Super High-Rise Building. Front Built Environ 1:. https://doi.org/10.3389/fbuil.2015.00011
  » https://doi.org/10.3389/fbuil.2015.00011
• Jacobsen E, Kootsookos P (2007) Fast, Accurate Frequency Estimators [DSP Tips & Tricks]. IEEE Signal Process Mag 24:123-125. https://doi.org/10.1109/MSP.2007.361611
  » https://doi.org/10.1109/MSP.2007.361611
• Jain VK, Collins WL, Davis DC (1979) High-Accuracy Analog Measurements via Interpolated FFT. IEEE Trans Instrum Meas 28:113-122. https://doi.org/10.1109/TIM.1979.4314779
  » https://doi.org/10.1109/TIM.1979.4314779
• Keliang D, Mingliang L, Zichao S, et al (2018) Application of Hilbert-Huang Transform in Health Inspection of Highway Bridge Structure. J Beijing Univ Civ Eng Archit 34:28-32
• Kunwar A, Jha R, Whelan M, Janoyan K (2013) Damage detection in an experimental bridge model using Hilbert-Huang transform of transient vibrations: DAMAGE DETECTION IN AN EXPERIMENTAL BRIDGE MODEL USING HHT. Struct Control Health Monit 20:1-15. https://doi.org/10.1002/stc.466
  » https://doi.org/10.1002/stc.466
• Mao Y-W, Tu Y-Q, Xiao W, Yang H-Y (2012) Advances and trends of study on discrete intensive frequency spectrum zooming analysis and correction methodology [J]. J Vib Shock 21:112-119+151
• Nikoo A, Roshandel Kahoo A, Hassanpour H, Saadatnia H (2016) Using a time-frequency distribution to identify buried channels in reflection seismic data. Digit Signal Process 54:54-63. https://doi.org/10.1016/j.dsp.2016.03.008
  » https://doi.org/10.1016/j.dsp.2016.03.008
• Rabiner L, Schafer R, Rader C (1969) The chirp z-transform algorithm. IEEE Trans Audio Electroacoustics 17:86-92
• Radad M, Gholami A, Siahkoohi HR (2015) S-transform with maximum energy concentration: Application to non-stationary seismic deconvolution. J Appl Geophys 118:155-166. https://doi.org/10.1016/j.jappgeo.2015.04.010
  » https://doi.org/10.1016/j.jappgeo.2015.04.010
• Sejdić E, Orović I, Stanković S (2018) Compressive sensing meets time-frequency: An overview of recent advances in time-frequency processing of sparse signals. Digit Signal Process 77:22-35. https://doi.org/10.1016/j.dsp.2017.07.016
  » https://doi.org/10.1016/j.dsp.2017.07.016
• Sharma M, Goyal D, Achuth PV, Acharya UR (2018) An accurate sleep stages classification system using a new class of optimally time-frequency localized three-band wavelet filter bank. Comput Biol Med 98:58-75. https://doi.org/10.1016/j.compbiomed.2018.04.025
  » https://doi.org/10.1016/j.compbiomed.2018.04.025
• Shen H, Li M, Luo F (2008) Strict maximum entropy spectral estimation based on recursive algorithm. Radar Sci Technol 4:288-291
• Wang Z, Han P, Cao J (2001) Overlapping Digital Filter Theory. SIGNAL Process. 17:189-191
• Xiao W, Dong W (2013) A frequency spectrum analysis method based on all-phase time-shift phase difference for themeasured signals of ship model test. J Ship Mech 17:998-1008
• Zhou J, Zhu Z, Shu Y, Cai Q (1999) Extrapolation of RF echo data based on AR modeling. Nanjing Univ Aeronaut Astronaut Trans 16:193-199

Edited by

Editor:

Marco L. Bittencourt.

Publication Dates

• Publication in this collection
  07 Aug 2020
• Date of issue
  2020

History

• Received
  01 May 2020
• Reviewed
  18 May 2020
• Accepted
  10 June 2020
• Published
  19 June 2020