Efficient Spectral Estimation of Non-Stationary Harmonic Signals Using Harmonic Transform

Non-stationary harmonic signals cannot be accurately represented by Fourier transform (FT). Fortunately, several methods for representing non-stationary harmonic signals exist including Fan-Chirp trasnform (FChT) or Harmonic transform (HT). This article is focused on the Harmonic transform and its computation. Estimation of slope of fundamental frequency change between analysed segments is essential for computation of the HT. The slope is estimated using several realisations of HT and comparing the spectral flatness. Optimisation of this procedure is presented in the article. Keywords—Harmonic transform, non-stationary harmonic signals, time-frequency representations.


I. INTRODUCTION
The FT is generally able to represent frequency content of a signal, when the signal is composed of components with invariant frequency.Such signals can be called stationary harmonic signals and by using FT we can get their frequency representation with sufficient resolution in a specified frequency band.The ability of the FT to represent frequency content of a signal dimnishes if the signals contains components with varying frequency [1], [2].
One solution of this problem is to use Warped Fourier Transform (WFT) [3], where the signal is frequency warped before applying the FT.This operation can be interpreted as change of the signal's scale for the conversion of time-varying frequency components to frequency invariant components.The scaling operation can be generalized using the Scale Transform (ST) [4], [5], [6], where the scale is taken as a physical property of the signal, or the scaling operation can be integrated into the definition of transformation, as in Harmonic Transform [7].
There are other transforms for representation of nonstationary harmonic signals, which are suitable for specific applications.The Fractional Fourier Transform (FrFT) [8], [9] and Chirp Transform (CT) [10], [11] are suitable for linearly changing frequency components, whereas the Fan-Chirp Transform [12], [13] is suitable for signals with frequency components varying linearly on a fan geometry.II.HARMONIC TRANSFORM Harmonic transform has been introduced in [7] and its main difference from Fourier transform is the integrated timewarping function.It is defined as M. Trzos and H. Khaddour are with the Department of Telecomunications, Brno University of Technology, Brno, Czech Republic, e-mail: trzos@phd.feec.vutbr.czand xkhadd00@stud.feec.vutbr.cz.
Manuscript received November 15, 2012; revised January 11, 2007.where φ u (t) is unit phase, which is the phase of the fundamental harmonic component divided by its instantaneous frequency [7], and φ u (t) is first derivation of φ u (t).Linear change of fundamental frequency in a given segment is presumed, which is sufficiently satisfied by selecting an analysis window of appropriate length.Instantaneous phase ϕ(t) of a sinusoid with linear change of frequency [14] is defined as where f 0 is fundamental fequency and = ∆f 0 /T is the change of fundamental frequency divided by lenght of the segment.Assuming discrete signal segment of the length N , where T = N/F s , the discrete phase ϕ(n) of a sinusoid with linear frequency variation [15] can be written as where f 0 is discrete instantaneous frequency, N is length of the analysis window, and F s is sampling frequency (initial phase is disregarded for simplicity).Initial fundamental frequency in a given segment can be written as where f c is the central fundamental frequency within a segment and a is the slope of fundamental frequency change within the segment.Substituting (4) to (3) we get [15] ϕ Frequencies of spectral lines of the Fourier transform are given as and from the equation ( 5) it is obvious that the instantaneous phase is Discrete harmonic transform (DHT) of signals with linear fundamental frequency variation [15] is defined as where is first-order derivation of (5).
The difference between DFT and DHT at analysis of non-stationary harmonic signals can be seen on a part of speech uttering from the PTDB-TUG database with frequency modulation.On Fig. 1 (bottom) we can see that the higher frequencies are smoothed due to frequency modulation.On Fig. 1 (top) even high frequency peaks are clearly visible.

III. ESTIMATING PARAMETER a
As can be seen from equation ( 4), a determines the change of fundamental frequency in the analysed signal, with the assumption of linear fundamental frequency change.Variation of a results in time-warping of the analysed signal, as can be seen on Fig. 2.
The parameter a is usually estimated by grid search.It consists of performing several realisations of HT with different a [7].Afterwards, the realisation which minimises spectral Comparison of spectral flatness measure and modified spectral measure fo the analysed segment.flatness measure (SFM) where DHT is discrete harmonic transform and |.| denotes absolute value, is selected as the correct value.This approach is computationaly extensive and is usually simplified restricting values a can have and by following the change of fundamental frequency in time [15].When using (10) the SFM has several minimums and if a search algorithm was used, it could fall into local minimum.It is also noteworthy that it is possible the harmonic transform |DHT(a, k)| will be equal to zero for some values of k, which could mean that the spectral flatness will be zero for all a. Removing zero values solves this problem and leads to bandlimited spectral flatness measure.
The harmonic spectrum is not complex conjugated even for real signals (which is true for Fourier transform).From the frequency axis point of view, the unit phase function φ u (t) shifts the spectrum towards lower frequencies if a is positive, and to higher frequencies if its negative.Using the formula ( 8) we get only one-sided spectrum, the right part will not represent harmonic components of the analysed signal.When estimating a, (10) has two minimums (see Fig. 3).For harmonic signal analysis, only left side of the spectrum is useful, because it appropriately represents non-stationary harmonic signal.Using the modified spectral flatness measure (MSFM) we can get function of a which has clearly defined minimum.This is caused by using only left side of the spectrum and it consequently leads to reducing the number of operations needed to compute spectral flatness by N 2 − 1.

IV. CONCLUSION
Harmonic transform with some of its features and parameters has been described.Furthermore, the parameter a and its influence on non-stationary signal analysis was analysed.By reducing the number of operations in computation of spectral flatness measure, the number of operations needed for computation of Harmonic transform is reduced.
Future work aims at further improving the performance of harmonic transform computation using parallelisation and using machine learning for estimation of parameter a. Sensitivity to noise will also be measured.
Fig.3.Comparison of spectral flatness measure and modified spectral measure fo the analysed segment.