A System-of-Systems Perspective on Frequency Estimation: Time-Frequency Distribution of Multiple LFM Signals

Abstract

This chapter provides a System-of-Systems (SoS) perspective on a study of frequency estimation of signals with a focus on Linear Frequency Modulation (LFM) signals. This chapter describes an SoS approach for frequency estimation using Chirplet Transform (CT), Hough Transform (HT), and the Short Time Fourier Transform (STFT) with filtering viewpoint. The filtering viewpoint employs the filter impulse response length to obtain the best time-frequency concentration for accurate estimation of a signal frequency. The optimum impulse response length can be found by varying the length of the filter impulse response and observe the changing in the time-frequency distribution (TFD). The chapter shows that when the length of the impulse response becomes longer, the time-frequency concentration in TFD increases first and then decreases.

1. Introduction

The study of stationary deterministic signal has been greatly explored and appreciated. On the other hand, most signals encountered in applications are random and nonstationary. Unlike the time-invariant statistical properties of stationary signal, the statistical properties of nonstationary signal are normally time-variant where a time-frequency combined analysis tool, time-frequency distribution (TFD), is required to observe the nonstationary signal in the time domain and the frequency domain at the same time. The early TFD is often given by Short Time Fourier Transform (STFT), Gabor Transform, or Continuous Wavelet Transform (CWT). The classical one is the quadratic Wigner-Ville Distribution (WVD). The latest type would be the parameterized TFD [1, 2, 3] developed in recent years. Both STFT and CWT, in the sense of transform, are not able to achieve a fine resolution in both time and frequency domain simultaneously, due to the restriction of the Heisenberg–Gabor inequality. Linear frequency modulated (LFM) is one of the pulse compression techniques in the Radar system to solve the conflict between rang and resolution where the carrier frequency is continuously modulated during the pulse duty time. The term of instantaneous frequency (IF) is used to describe how the carrier start frequency changes linearly all the way up to the end frequency. In fact, the LFM signal has no fixed period nor frequency within each pulse duty time. The quadratic WVD will achieve the highly accurate frequency component for noise-free LFM signal, where the constant amplitude brings WVD a row of delta functions along the linear IF trajectory [4, 5]. In the case of noisy LMF, the WVD peak position will bias from the true IF, where the bias-to-variance tradeoff is inevitable in the IF estimation. The Chirplet Transform (CT) [1, 2] is a typical parametric TFD, which is particularly designed for the analysis of chirp-like signals with linear IF. In the initialization process of the CT, the parameters estimation is based on the peak of the STFT magnitude, thus the estimation results are greatly affected by the background noise. For the multiple LFM signal, it is difficult to distinguish and track multiple IF lines. If the Hough Transform (HT) is applied to the spectrogram magnitude first, the robust parameter estimation can be obtained for each component, a set of time-frequency images can then be emerged by post processing to finally get a TFD with higher concentration.

Some scholars have analyzed the nonstationary signal with filtering viewpoint. In order to analyze the audio signal, Brown proposed the constant Q transform (CQT) [6], where the central frequency of each band is not uniformly distributed and its frequency resolution is not a fixed value in the frequency domain, that is more suitable for nonstationary audio signal processing. Another adaptive filter bank is proposed in [7] where the frequency resolution is changed by adjusting the window length in each sub-band. A more generalized TFD is proposes based on the traditional CQT [8]. It can be used to define a time-frequency analysis framework with arbitrary central frequency at arbitrary frequency resolution. The parameters in the framework are clearly defined to achieve a good resolution at any given frequency range. Another novel time-frequency analysis is proposed in [9] where the filter bank is a high-resolution Gaussian filter bank. Based on the nonlinear characteristics of the human auditory system, the Gaussian filter bank is designed to adjust the central frequency of each band. At the same time, the multi-resolution characteristic of the filter bank is discussed based on the idea of Wavelet transform.

However, little has been reported about the STFT in the filtering viewpoint. The STFT is always regarded as time shifted Fourier transformations where frame length is fixed for every transformation. If STFT is treated as the outputs of a filter bank, the length of each channel impulse response can be set differently according to the different signal frequency of different channel. This chapter combines our two recent conference papers [10, 11] and provides a unified framework derived from a System-of-Systems perspective for analyzing LFM signal using Filtering Viewpoint approach. Based on a series of experiments, the impact of the filter impulse response length is observed with the variation of TFD. It is also proved that longer filter does not always guarantee a better energy concentration in the TFD. To obtain the best TFD for the LFM signal processing, an optimal impulse response length needs to be determined beforehand.

This chapter presents a SoS approach for frequency estimation with a focus on LFM signal. Using a standard system engineering approach, we can decompose the frequency estimation process into three systems consisting of:

• System 1: responsible for transforming the time-domain signal into (i) frequency-domain signal, and (ii) TF transform using CWT. For LFM signal type, CT is selected for CWT. Section 2 provides detailed description of STFT and CT transforms.

• System 2: responsible for detecting instantaneous signal frequency (IF) in the presence of noise. Section 3 describes a proposed technique using Hough transform for detecting IF straight lines.

• System 3: responsible for (i) assessing of filter impulse response length on TFD (see Section 4), and (ii) analyzing the time-frequency behavior and applying Hough transform for frequency estimation (see Section 5).

Section 4 presents a series of experiments showing the influence of the impulse response length to the time-frequency concentration in TFD and provides the steps of finding the optimal impulse response length. Section 5 provides an example of time-frequency analysis and proves the feasibilities of the proposed approach before the conclusion.

2. Frequency transforms: a review

2.1 STFT using filter bank viewpoint

The STFT is normally regarded as the Fourier transform of the framed signals with an observation window of fixed length. Whereas being looked from the filtering viewpoint, the filter bank has some advantages that the traditional transform does not have. For example, if the frequency range of the signal is known, the corresponding bands can be selected in advance and only the selected bands need to be calculated, which will greatly reduce the computation cost. From the viewpoint of filtering, STFT can be actually regarded as the outputs of a filter bank, in which each band is called “analytical filter” and has its own impulse response function, w [n]. The STFT at the band centered at ω₀ is the output of the analytical filter driven by the signal demodulated with the complex exponential carrier at the central frequency of that band.

$$X\left(n,{\omega }_0\right)=\sum _{m=-\infty }^{\infty }\left(x\left[m\right]e^{-j{\omega }_{0}m}\right)w\left[n-m\right]$$E1

Rewrite the demodulated signal into x_w [n], the STFT at ω₀ has no difference to the discrete convolution between the driven signal x_w[n] and impulse response w[n].

$$X\left(n,{\omega }_0\right)=x_w\left[n\right]\ast w\left[n\right]$$E2

This can also be expressed as:

$$X\left(n,{\omega }_0\right)=e^{-j{\omega }_{0}n}\left(x\left[n\right]\ast w\left[n\right]e^{j{\omega }_{0}n}\right)$$E3

where the sequence x[n] passes a filter with the modulated impulse response first, before getting demodulated with the same carrier.

Following this idea, the entire STFT covering N discrete DFT frequencies can be obtained as the output of the filter bank shown in Figure 1, where K = N-1, and N is the frame length of STFT. In fact, when this system deployed directly into the hardware, the carrier frequency of each band-pass filter can be arbitrarily selected, which could have nothing to do with N. Another benefit of the direct deployment is the time resolution of this filter bank output, which could be as high as the time resolution of the input signal x[n]. The challenge is the selection of response time of each selected band may be different from band to band because of the different filter lengths are required to provide high time-frequency concentration in the filter bank output x[n, k].

Figure 1.

STFT at N DFT frequencies using filtering viewpoint.

2.2 Chirplet transform (CT)

The CT of a signal s(t) is defined as

$${\mathit{CT}}_s\left(t_0,\omega ,\alpha ;\sigma \right)=\underset{-\infty }{\overset{\infty }{\int }}z\left(t\right){\psi }_{t_0,\alpha ,\sigma }^{\ast }\left(t\right)exp\left(-\mathit{j\omega t}\right)\mathit{dt}$$E4

where z(t) is the analytical signal of s(t) generated by the Hilbert transform, and ψ is a complex window given by

$${\psi }_{t_0,\alpha ,\sigma }\left(t\right)=w_{\sigma }\left(t-t_0\right)exp\left[-j\frac{\alpha }{2}{\left(t-t_0\right)}^2\right]$$E5

where t₀ is the time shift and α the chirp rate. The window function w is usually taken as Gaussian function expressed as

$$w_{\sigma }\left(t\right)=\frac{1}{\sigma \sqrt{2\pi }}exp\left[-\frac{1}{2}{\left(\frac{t}{\sigma }\right)}^2\right]$$E6

From this definition, it can be seen that the CT can be decomposed into a series of operations: 1) rotating the signal under consideration by an angle in the time-frequency plane; 2) shifting the signal by a frequency increment; and 3) applying STFT with the Gaussian window.

This process can be depicted in the Figure 2. The solid line is the IF line of the target LFM signal that has the IF function ω(t) = ω₀ + λ₀t. The dot-dashed line is the IF line after the rotation, and the dashed line represents the IF line after translation.

Figure 2.

The Chirplet transform.

As mentioned earlier, given a set of properly determined kernel characteristic parameters, the CT could produce a high-quality TFD for a considered signal. The result can have an excellent T-F concentration, which measure the IF trajectory width over the TFD surface, so the IF trajectory can be easily identified. Therefore, the determination of proper parameters is critical for the application of the CT method. Briefly speaking, the basic idea of the CT based T-F analysis uses the kernel characteristic parameter (α = 0) to form the TFD, and then finds the maximum value along time axis in the time-frequency plane. The resulting maximum line approximation is considered to be an IF trajectory. The chirp parameter obtained by the line fitting will be reapplied to the CT transformation. The procedure can be repeated until no evident improvement is observed in the T-F concentration of the TFD.

To measure the T-F concentration of the TFD, the Rényi entropy can be used with the definition:

$$E\left(s\right)=-\int \int log{\left|{\mathit{CT}}_s\left(t,\omega \right)\right|}^3\mathit{\text{dtdω}}$$E7

And the termination condition can be set as

$${\xi }_s=\left|E_{i+1}\left(s\right)-E_i\left(s\right)\right|<\delta$$E8

or

$${\xi }_s=\frac{\left|E_{i+1}\left(s\right)-E_i\left(s\right)\right|}{\left|E_{i+1}\left(s\right)\right|}<\delta$$E9

where δ is a predetermined threshold. If the parameters of the initialization are not accurate, it will lead to be a lower overall T-F concentration. Because STFT is susceptible to noise, the robustness of the whole algorithm is not strong. Fortunately, the fractional Fourier transform can compensate for this shortcoming.

3. Detect IF straight lines by Hough transform

Through the description of the previous section, we know that as long as we can get the precise FM parameters, we can get the TFD with high T-F concentration. However, for multicomponent signal with low SNR, the difficulty will be significantly huge. Inspired by the reference [12], we apply the robust HT to detect the IF lines in the CT-TFD to depress the noise during the process of line fitting. The reason why we adopt this image processing technique is that the HT can detect multiple lines accurately, even in the low SNR situations.

For an IF straight line in the Cartesian coordinate plane, there are two common representations: point-slope form and two points form. In the HT, however, another representation is considered: coordinate (r, θ) is used to represent a straight line where r is the distance from the straight line to the origin, θ is the angle between the x-axis and the perpendicular line passing through the origin as shown in Figure 3.

Figure 3.

The line representation in Hough transform.

The idea of using HT to detect straight lines needs to assume M straight lines for each point, usually M = 180. In this chapter, the angle resolution of the detected line is 1 degree. The coordinates (r, θ) of the M lines are calculated respectively. If there are a total of L points to be checked, the total numbers of coordinates are M × L. If a number of points are on the certain same line, then there must be a same number of points with the same the same coordinate (r_c, θ_c). The following example shows that if there are three points in the plane, we can determine whether or not the three points are on the same straight line by their HT coordinates.

4. Impact of filter impulse response length on TFD: simulation results

The below simulations are all conducted on the Matlab 2015b installed in Windows10 system on the Dell T7910 workstation with 2 Intel Xeon E5-2630v3 CPUs and 256G LRDIMM memories.

To see the influence of the filter impulse response length on the concentration of the TFD, the LFM signals of different FM parameters start with the simulated signal:

$$z_1\left(t\right)=exp\left[j\left(12\mathit{\pi t}+10\pi \right)t\right]$$E10

where its IF line function is f₁(t) = 6 t + 5. The signal lasts for 5 seconds, and 1000 channels are selected uniformly from 0 ∼ 50 Hz. The sampling frequency of z₁(t) is fixed at 100 Hz, but the impulse response length of all the bands is tested on 40, 80, 100 and 150 samples respectively. The corresponding TFDs are shown by Figure 4.

Figure 4.

TFDs of z₁(t) using filter bank at 4 different filter lengths of 40 (a), 80 (b), 100(c), and 150 (d) samples.

It can be seen clearly that as the impulse response length grows, the T-F concentration of the TFD increases first and then decreases. In order to check whether this is a universal phenomenon, the following multiple LFM signal is built.

$$z_2\left(t\right)=\sum _{k=0}^{7}exp\left[j\left(2k+1\right)5\pi t^2+10\mathit{\pi t}\right]$$E11

The IF slopes of all the components of z₂(t) form a sequence of common difference of 5 Hz, starts from 2.5 Hz to the highest 37.5 Hz. Thus, the IF of z₂(t) is in the range of 5 ∼ 192.5 Hz within the 5 s duration. In the Figure 5, the TFDs are given for the impulse response length of 50, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550 and 600 samples respectively at the sampling frequency of 1000 Hz. 300 bands are uniformly selected from 0 ∼ 300 Hz in the calculation of each TFD. It is shown clearly again by the Figure 5, as the impulse response length increases, the T-F concentration of each component in the signal z₂(t) increases first and then decreases. The larger the frequency modulation coefficient (or IF slope) is, the faster its TFD reaches its best time-frequency concentration. For the component with the lowest IF slope of 2.5 Hz/s, the concentration change is not obvious since it is close to a stationary signal.

Figure 5.

TFDs of z₂(t) using filter bank with 12 different filter lengths running from 50 (a) to 600 (l) samples at the step of 50 samples.

In order to verify the above conclusion, a complex stationary signal below of 5 Hz is also tested.

$$z_3\left(t\right)=exp\left(j10\mathit{\pi t}\right)$$E12

The testing range of the impulse response length, the selected channels, the signal duration, and the sampling frequency are taken in the same way as that of z₁(t). The resulted TFDs are shown in the Figure 6. It can be seen that for this stationary signal z₃(t), as the impulse response length increases, the T-F concentration of the TFD keeps increasing because longer filter collects nothing but more energy of the stable signal, which is different from the case of LFM signal where more interferences will be observed by the longer filter. This is to say that the length of the impulse response should match the changing speed of the signal IF, so that the TFD could accurately reflect the T-F trajectories of the LFM signals.

Figure 6.

TFDs of z₃(t) using filter bank at 4 different filter lengths of 40 (a), 80 (b), 100(c), and 150 (d) samples.

The next simulation is about the harmonically related multiple LFM signals as given below.

$$z_4\left(t\right)=\sum _{k=0}^{5}exp\left[j6\pi t^2+20\left(2k+1\right)\mathit{\pi t}\right]$$E13

The length of the impulse response runs within {150, 250, 350, 450, 550, 650}. The selected channels, the signal duration, and the sampling frequency are taken in the same way as what has been done on z₂(t). The resulted TFDs are shown in the Figure 7.

Figure 7.

TFD magnitudes of z₄(t) using filter bank at the filter length of 150 (a), 250 (b), 350 (c), 450 (d), 550 (e), and 650 (f) samples.

It can be seen from Figure 7 that with the increasing of impulse response length, the TFD T-F concentration of each component of the signal z₄(t) gets better and better, while the T-F concentrations of different components with the same filter length have no difference from each other. In other words, the same IF slope results the same TFD T-F concentration given the same filter length.

Through the description in the previous section, the high T-F concentration over the TFD surface can be obtained if the impulse response length matches to the IF slope of FLM signal. In order to obtain this optimal filter length, the T-F concentration M(l) is measured by the averaged bandwidth in the TFD as given below, where A_l(t,ω) is the TFD magnitude calculated with the filter length of l at the (t,ω) point, T is the time limit of TFD. Based on this TFD concentration measurement, the optimal filter length can be determined by Eq. (15). The IF of the signal of interests can then be accurately estimated from the resulted A_Lopt(t,ω) using any modest ridge detection method.

$$\begin{array}{l}M\left(l\right)=\frac{1}{T}\sum _{t=1}^T\sum _{\omega }{B}_l\left(t,\omega \right)\\ B_l\left(t,\omega \right)=\left\{\begin{array}{cc}1,& A_l\left(t,\omega \right)>0.9\underset{arg\omega }{max}\left[A_l\left(t,\omega \right)\right]\\ 0,& \mathit{\text{otherwise}}\end{array}\right.\end{array}$$E14

$$L_{\mathit{opt}}=\underset{l}{argmax}\left[M\left(l\right)\right]$$E15

When the LFM slope coefficient is given, the channels of the filter bank can be precisely selected to cover that IF range. Even for the signal with unknown LFM slope coefficient, one can always observe the rough IF range using traditional STFT, so no difficulty will be met in the filter bank channel selection. Figure 8 compares the proposed TFD of filtering viewpoint with the traditional STFT of transform viewpoint. The LFM signal is z₁(t) sampled at the frequency of 100 Hz, the impulse response length is 80 samples for all the channels, the selected 1000 channels uniformly divide the range of 0 ∼ 50 Hz, which calls for the frame length of 2000 samples in the traditional STFT to reach the same frequency resolution. It can be seen that the concentration of the proposed TFD is significantly better than that of the traditional STFT.

Figure 8.

Proposed TFD (a) vs. traditional STFT (b) of z₁(t) at the same frequency resolution.

5. Frequency estimation: time-frequency analysis for multiple LFM signal based on the CT and HT

Taking the two-component signal as an example, it will give more accurate parameters through line fitting before applying them to CT. In order to see the TFD difference between parameter matched and nonmatched CTs, a simulated multiple LFM signal is considered as:

$$s\left(t\right)=sin\left[2\pi \left(40t-t^2\right)\right]+sin\left[2\pi \left(10t+\frac{5}{4}{t}^2\right)\right]\left(0\le t\le 15s\right)$$E16

where IF lines are f₁(t) = 40-2 t and f₂(t) = 10 + 2.5 t. So, the actual parameters should be α₁ = −4π and α₂ = 5π. Suppose we have these two exact parameters, and then apply them into the CT transformation and observe the properties of the TFD graphs.

(a) α₁ = −4π, (b) α₂ = 5π, (c) the superposition of (a) and (b).

From Figure 9, one can see that when the parameter α = α₁, the component with IF of f₁(t) has better T-F concentration in the CT-TFD while the component with IF of f₂(t) has very low T-F concentration and also very low magnitude. It shows the opposite situation when α = α₂. Adding the two spectra together, a good TFD can be obtained as shown above. If we cut off the very low magnitude with a threshold before adding up the two spectra, a sharper TFD with higher T-F concentration can be obtained as shown by Figure 10.

Figure 9.

The TFD results of s(t) using CT with parameter (a) α₁ = -4π, (b) α₂ = 5π, (c) the superposition of (a) and (b).

Figure 10.

The results of threshold filtering: (a) α₁ = -4π, (b) α₂ = 5π, (c) the superposition of (a) and (b).

In order to highlight the advantages of HT, the analog signal of Eq. (16) is generated at the SNR of 2 dB, 0 dB, and -2 dB. Then the STFT and the corresponding HT are calculated as shown in Figure 11. From the below results, one can see clearly that in the case of low SNR, the STFT is relatively fuzzy, and the ridge edge extraction alone will not give the good results. However, the accuracy of HT is relatively high at all the 3 SNR levels. That is the reason why the HT is adopted in the time-frequency analysis of multiple LFM signal.

Figure 11.

The STFT (a), the corresponding HT (b), and the TFD after CT & HT (c) of s(t) at the SNR of 2 dB(top), 0 dB(middle), and -2 dB(bottom).

6. Conclusions

This chapter investigates an SoS approach for frequency estimation using TFD calculation techniques of the LFM signal through a linear filtering viewpoint. The influence of filter length on the TFD concentration is closely observed through a series of simulations. The simulation results show that the IF slope of the LFM signal is related with the optimal filter length, the higher the slope is, the shorter the optimal filter length is. On the other hand, the same IF slope results the same TFD concentration at the same filter length, no matter how high the IF is. Under the same time-frequency resolution, the traditional STFT shows significantly lower time-frequency concentration than that of the TFD obtained by the proposed filter bank based on the filtering viewpoint. Thought the channel frequency can be freely selected as needed, which saves the computation in the irrelevant frequency bands, many nonstationary signals are not always linearly modulated. For the signals with the nonlinear FM coefficient, the advantage of the filter bank TFD is no longer obvious. This is also a question that needs to be studied further.

The CT has the advantage of high T-F centralization but is easily affected by the noise. In addition, many nonstationary signals are of multicomponent. The decomposition of multicomponent signal [12, 13] into single component signals under noise conditions is a difficult problem. For the multicomponent LFM signal, the Hough transform is adopted to the parametric T-F analysis to obtain the result with good concentration. However, for the multicomponent nonlinear FM signal, there is no effective trajectory detection method, so it will be more difficult to decompose each nonlinear FM component. That will be the focus of our future work.

Acknowledgments

The authors would like to express great thanks to the Shanghai Key Lab in Information Security Management (AGK201709), and Shandong Nature Science Fund Committee (ZR2016FM44) for their financial supports.