Any wave of complex shape can be represented as the sum of simple waves.
Joseph Fourier was keen to describe in mathematical terms how heat travels through solid objects ( cm. Heat exchange). Perhaps his interest in warmth flared up while he was in North Africa: Fourier accompanied Napoleon on a French expedition to Egypt and lived there for some time. To achieve his goal, Fourier had to develop new mathematical methods. The results of his research were published in 1822 in the work "Analytical Theory of Heat" ( Theorie analytique de la chaleur), where he told how to analyze complex physical problems by decomposing them into a number of simpler ones.
The method of analysis was based on the so-called Fourier series. In accordance with the principle of interference, the series begins with the decomposition of a complex shape into simple ones - for example, a change in the earth's surface is due to an earthquake, changes in the orbit of a comet due to the influence of the attraction of several planets, a change in heat flow due to its passage through an irregularly shaped obstacle made of heat-insulating material. Fourier showed that a complex waveform can be represented as the sum of simple waves. As a rule, the equations describing classical systems are easily solved for each of these simple waves. Fourier went on to show how these simple solutions can be summarized to get the solution of the entire complex problem as a whole. (Mathematically speaking, a Fourier series is a method of representing a function as a sum of harmonics—sine and cosine, so Fourier analysis was also known as harmonic analysis.)
Until the advent of computers in the mid-twentieth century, Fourier methods and the like were the best weapons in the scientific arsenal when attacking the complexities of nature. Since the advent integrated methods Fourier scientists were able to use them to solve not only simple problems that can be solved by direct application of Newton's laws of mechanics and other fundamental equations. Many of the great achievements of Newtonian science in the 19th century would in fact have been impossible without the use of the methods first proposed by Fourier. Subsequently, these methods were applied in solving problems in various fields from astronomy to mechanical engineering.
Jean-Baptiste Joseph Fourier
Jean-Baptiste Joseph Fourier, 1768-1830
French mathematician. Born in Auxerre; at the age of nine he was left an orphan. Already at a young age he showed aptitude for mathematics. Fourier was educated at a church school and a military school, then worked as a teacher of mathematics. Throughout his life he was actively involved in politics; was arrested in 1794 for defending the victims of terror. After the death of Robespierre, he was released from prison; took part in the creation of the famous Polytechnic School (Ecole Polytechnique) in Paris; his position provided him with a springboard to advance under Napoleon's regime. Accompanied Napoleon to Egypt, was appointed governor of Lower Egypt. Upon returning to France in 1801, he was appointed governor of one of the provinces. In 1822 he became permanent secretary of the French Academy of Sciences, an influential position in the scientific world of France.
FOURIER TRANSFORM AND CLASSICAL DIGITAL SPECTRAL ANALYSIS.
Medvedev S.Yu., Ph.D.
Introduction
Spectral analysis is one of the signal processing methods that allows one to characterize the frequency composition of the measured signal. The Fourier transform is a mathematical framework that relates a temporal or spatial signal (or some model of this signal) to its representation in the frequency domain. Statistical methods play an important role in spectral analysis, since signals are usually random or noisy during propagation or measurement. If the main statistical characteristics of a signal were precisely known, or could be determined from the finite interval of this signal, then spectral analysis would be a branch of "exact science". However, in reality, only an estimate of its spectrum can be obtained from a signal segment. Therefore, the practice of spectral analysis is a kind of craft (or art?) of a rather subjective nature. The difference between the spectral estimates obtained as a result of processing the same signal segment by different methods can be explained by the difference in the assumptions made regarding the data, different methods of averaging, etc. If the characteristics of the signal are not known a priori, it is impossible to say which of the estimates is better.
Fourier transform - the mathematical basis of spectral analysis
Let us briefly discuss the different types of Fourier transform (for more details, see ).
Let's start with the Fourier transform of a time-continuous signal
,
(1)
which identifies the frequencies and amplitudes of those complex sinusoids (exponential), into which some arbitrary oscillation is decomposed.
Reverse transformation
. (2)
The existence of the forward and inverse Fourier transform (which we will call the continuous-time Fourier transform - CTFT) is determined by a number of conditions. Sufficient - absolute signal integrability
. (3)
A less restrictive sufficient condition is the finiteness of the signal energy
. (4)
We present a number of basic properties of the Fourier transform and the functions used below, noting that the rectangular window is defined by the expression
(5)
and the sinc function is an expression
(6)
The function of samples in the time domain is determined by the expression
(7)
This function is sometimes also called the periodic continuation function.
Table 1. The main properties of the NVPF and functions
|
Property, function |
Function |
transformation |
| Linearity |
ag(t) + bh(t) |
aG(f) + bH(f) |
| Timeshift |
h (t - t0) |
H(f)exp(-j2pf t 0) |
| Frequency shift (modulation) |
h (t)exp(j2pf0 t) |
H(f - f0) |
| Scaling |
(1 / |a|)h(t / a) |
H(af) |
| Time domain convolution theorem |
g(t)*h(t) |
|
| Convolution theorem in the frequency domain |
g(t) h(t) |
G(f)*H(f) |
| window function |
Aw(t / T) |
2ATsinc(2Tf) |
| sinc function |
2AFsinc(2Ft) |
Aw(f / F) |
| impulse function |
Ad(t) |
|
| Counts function |
T(f) |
FF(f), F=1/T |
Another important property is established by Parseval's theorem for two functions g(t) and h(t):
. (8)
If we put g(t) = h(t), then Parseval's theorem reduces to the theorem for the energy
. (9)
Expression (9) is, in essence, just a formulation of the energy conservation law in two domains (time and frequency). In (9), the total energy of the signal is on the left, so the function
(10)
describes the frequency distribution of energy for a deterministic signal h(t) and is therefore called the energy spectral density (SPD). Expressions
(11)
one can calculate the amplitude and phase spectra of the signal h(t).
Discretization and weighting operations
In the next section, we will introduce the Discrete Time Fourier Series (DTFT) or otherwise the Discrete Fourier Transform (DFT) as a special case of the Continuous Time Fourier Transform (CTFT) using two basic signal processing operations - sampling ( discretization) And weighing using a window. Here we consider the influence of these operations on the signal and its transformation. Table 2 lists the functions that perform weighting and discretization.
With uniform readings with an interval of T seconds, the sampling frequency F is equal to 1 /T Hz. Note that the weighting function and sampling function in the time domain are denoted respectively by TW (time windowing) and TS (time sampling), and in the frequency domain by FW (frequency windowing) and FS (frequency sampling).
Table 2. Weighting and discretizing functions
|
Operation |
Time function |
transformation |
| Time domain weighting (window width NT sec) |
TW=w(2t / NT - 1) |
F(TW)=NTsinc(NTf)exp(-jpNTf) |
| Frequency domain weighting (window width 1/T Hz) |
|
FW=w(2Tf) |
| Time readings (interval T sec) |
TS=T T (t) |
|
| Frequency sampling (1/NT Hz interval) |
|
|
Let us assume that samples are taken of a continuous real signal x(t) with a limited spectrum, the upper frequency of which is equal to F0. The NITF of a real signal is always a symmetrical function with a total width of 2F0, see Fig.1.
Signal samples x(t) can be obtained by multiplying this signal by the function of samples:
(12)
Fig.1 is an illustration of the sampling theorem in the time domain for a real spectrum limited signal:
a - the original function of time and its Fourier transform;
b - the function of readings in time and its Fourier transform;
c - time readings of the original function and its periodically extended Fourier transform for the case Fo<1/2T;
d - frequency window (ideal filter low frequencies) and its Fourier transform (sinc function);
d is the original time function, restored by means of a convolution operation with the sinc function.
According to the frequency domain convolution theorem, the CTF of the signal x(t) is simply the convolution of the spectrum of the signal x(t) and the Fourier transform of the time sample (TS) function:
. (13)
The convolution X(f) with the Fourier transform of the sampling function F (TS)=Y1/T(f) simply continues X(f) periodically with a frequency interval of 1/T Hz. Therefore, XS(f) is a periodically extended spectrum of X(f). In the general case, samples in one domain (for example, time) lead to a periodic continuation in the transformation domain (for example, frequency). If the sampling frequency is chosen low enough (F< 2Fo), то периодически продолженные
спектры будут перекрываться с соседними. Это перекрытие носит название
эффекта наложения в частотной области.
In order to restore the original time signal from its readings, i.e. to interpolate some continuum of values between these samples, it is possible to pass the sampled data through an ideal low-pass filter with a rectangular frequency response (Fig. 1d)
. (14)
As a result (see Fig. 1 e) the original Fourier transform is restored. Using the convolution theorems in time and frequency domains, we get
. (15)
Expression (15) is a mathematical notation sampling theorems in the time domain(theorems of Whittaker, Kotelnikov, Shannon - UKSH), which states that with the help of the interpolation formula (15), a real signal with a limited spectrum can be exactly restored by an infinite number known time readings taken with frequency F і 2F0. The dual to theorem (15) is the theorem samples in the frequency domain for signals with a limited duration.
Operations in the time domain similar to (14) are described by the expression
, (16)
and the corresponding transformations - by expressions
Thus, the NITF X(f) of a certain signal with a limited duration can be unambiguously restored from equidistant samples of the spectrum of such a signal, if the selected frequency sample interval satisfies the condition F1/2T 0 Hz, where T 0 is the signal duration.
Relations between continuous and discrete transformations
A pair of transformations for the usual definition of the Discrete Fourier Transform (DFT) of an N-point time sequence x[n] and the corresponding N-point Fourier transform sequences X[k] is given by the expressions
,
(18)
. (19)
In order to obtain spectral estimates from data samples in the corresponding units of energy or power, we write the discrete-time Fourier series (DTFT), which can be considered as some approximation of the continuous-time Fourier transform (CTFT), based on the use of a finite number of data samples:
In order to show the nature of the correspondence to the dvrf ( discrete functions in both the time and frequency domains) and CTF (continuous functions in the time and frequency domains), we need a sequence of four linear commutation operations: weighting in the time and frequency domains, and sampling or sampling both in the time and frequency domains. If the weighting operation is performed in one of these areas, then, according to the convolution theorem, it will correspond to the execution of the filtering (convolution) operation in another area with the sinc function. Similarly, if discretization is performed in one region, then the periodic continuation operation is performed in the other. Since weighing and sampling are linear and commutative operations, there are various ways of ordering them, giving the same final result with different intermediate results. Figure 2 shows two possible sequences for these four operations.
Rice. 2. Two possible sequences of two weighing operations and two sampling operations linking the NIPF and DTRF: FW - application of a window in the frequency domain; TW - application of a window in the time domain; FS - sampling in the frequency domain; TS - sampling in the time domain.
1 - Fourier transform with continuous time, equation (1);
4 - Fourier transform with discrete time, equation (22);
5 - Fourier series with continuous time, equation (25);
8 - Fourier series with discrete time, equation (27)
As a result of performing the weighing and sampling operations at nodes 1, 4, 5 and 8, there will be four different types of Fourier relations. Nodes where the function in frequency domain is continuous, refer to transformations Fourier, and the nodes in which the function is in the frequency domain discrete refer to Fourier series(for more details, see).
So at node 4, weighting in the frequency domain and sampling in the time domain generates discrete-time transformation Fourier transform (DTFT), which is characterized by a periodic function of the spectrum in the frequency domain with a period of 1/T Hz:
(22)
(23)
Note that expression (22) defines a certain periodic function that coincides with the original transformed function specified in node 1 only in the frequency range from -1/2T to 1/2T Hz. Expression (22) is related to the Z-transform of the discrete sequence x[n] by the relation
(24)
So the DTFT is just the Z-transform computed on the unit circle and multiplied by T.
If we move from node 1 to node 8 in Fig. 2 along the lower branch, at node 5, the operations of weighting in the time domain (limiting the signal duration) and sampling in the frequency domain generate a continuous-time Fourier series (CTSF). Using the properties and definitions of functions given in tables 1 and 2, we obtain the following pair of transformations
(25)
(26)
Note that expression (26) defines a certain periodic function that coincides with the original one (at node 1) only on the time interval from 0 to NT.
Regardless of which of the two sequences of four operations is chosen, the final result at node 8 will be the same - discrete-time Fourier series, which corresponds to the following pair of transformations obtained using the properties indicated in Table 1.
, (27)
where k=-N/2, . . . ,N/2-1
,
(28)
where n=0, . . . ,N-1 ,
The energy theorem for this WWRF is:
,
(29)
and characterizes the energy of a sequence of N data samples. Both sequences x[n] and X[k] are periodic modulo N, so (28) can be written in the form
, (30)
where 0 n N. The factor T in (27) - (30) is necessary so that (27) and (28) are actually an approximation of the integral transformation in the integration domain
.(31)
Zero padding
Through a process called zero padding, the discrete time Fourier series can be modified to interpolate between N values of the original transform. Let the available data samples x,...,x be supplemented with zero values x[N],...X. The TDRF of this zero-padded 2N-point data sequence will be given by
(32)
where the upper sum limit on the right is modified to accommodate null data. Let k=2m, so
,
(33)
where m=0,1,...,N-1, defines even values of X[k]. This shows that for even values of the index k, the 2N-point discrete-time Fourier series reduces to an N-point discrete-time series. Odd values of index k correspond to interpolated TDGF values located between the values of the original N-point TDGF. As more and more zeros are added to the original N-point sequence, even more interpolated data can be obtained. In the limiting case of an infinite number of introduced zeros, the DTRF can be considered as a discrete-time Fourier transform of an N-point data sequence:
.
(34)
Transformation (34) corresponds to node 6 in Fig.2.
There is a misconception that zero-padding improves resolution because it increases the length of the data sequence. However, as follows from Fig. 3, padding with zeros does not improve the resolution of the transformation obtained from the given finite data sequence. Zero padding simply allows you to get an interpolated transformation more flattened shape. In addition, it eliminates the uncertainties due to the presence of narrowband signal components whose frequencies lie between N points corresponding to the estimated frequencies of the original TPDF. When padded with zeros, the accuracy of estimating the frequency of spectral peaks also increases. By the term spectral resolution we mean the ability to distinguish between the spectral responses of two harmonic signals. A generally accepted rule of thumb, often used in spectral analysis, states that the frequency separation of distinguishable sinusoids cannot be less than equivalent window bandwidth, through which segments (segments) of these sinusoids are observed.
Fig.3. Interpolation by padding with zeros:
a - DPRF module for 16-point data recording containing three sinusoids without zero padding (uncertainties are visible: it is impossible to say how many sinusoids are in the signal - two, three or four);
b - TDWF module of the same sequence after a twofold increase in the number of its readings due to the addition of 16 zeros (uncertainties are resolved, since all three sinusoids are distinguishable;
c - the TDWF module of the same sequence after a fourfold increase in the number of its readings due to the addition of zeros.
The equivalent window bandwidth can be defined as 
where W(f) is the discrete-time Fourier transform of the window function, for example, rectangular (5). Similarly, you can enter equivalent window duration
It can be shown that the equivalent duration of a window (or any other signal) and the equivalent bandwidth of its transformation are mutually reciprocal: TeBe=1.
Fast Fourier Transform
The Fast Fourier Transform (FFT) is not just another variation of the Fourier transform, but rather the name of a range of efficient algorithms, designed for fast calculation of the discrete-time Fourier series. The main problem that arises in the practical implementation of the WWRF lies in the large number of computational operations proportional to N2. Although long before the advent of computers, several efficient computational schemes were proposed that could significantly reduce the number of computational operations, a real revolution was made by the publication in 1965 of an article by Cooley (Cooly) and Tukey (Tukey) with a practical algorithm for fast (number of operations Nlog 2 N) calculation of the DTWF . After that, many variants, improvements and additions to the basic idea were developed, which made up the class of algorithms known as the fast Fourier transform. The basic idea behind the FFT is to divide an N-point TDGF into two or more TDGFs of smaller length, each of which can be computed separately and then linearly summed with the others to obtain the TDFT of the original N-point sequence.
We represent the discrete Fourier transform (DTFT) in the form
, (35)
where the value of W N =exp(-j2 /N) is called the turning factor (hereinafter in this section, the sampling period is T=1). Select from the sequence x[n] elements with even and odd numbers
. (36)
But since that
. Therefore, (36) can be written as
, (37)
where each of the terms is a transformation of length N/2
(38)
Note that the sequence (WN/2) nk is periodic in k with period N/2. Therefore, although the number k in expression (37) takes values from 0 to N-1, each of the sums is calculated for values of k from 0 to N/2-1. One can estimate the number of complex multiplication and addition operations required to calculate the Fourier transform in accordance with the algorithm (37)-(38). Two N/2-point Fourier transforms according to formulas (38) require 2(N/2) 2 multiplications and approximately the same number of additions. The union of two N/2-point transformations according to formula (37) requires N more multiplications and N additions. Therefore, to calculate the Fourier transform for all N values of k, it is necessary to perform N+N 2 /2 multiplications and additions. At the same time, direct calculation by formula (35) requires N 2 multiplications and additions. Already for N>2, the inequality N+N 2 /2< N 2 , и, таким образом, вычисления по алгоритму (37)-(38) требуют меньшего числа математических операций по сравнению с прямым вычислением преобразования Фурье по формуле (35). Так как вычисление N-точечного преобразования Фурье через два N/2-точечных приводит к экономии вычислительных операций, то каждое из N/2-точечных ДПФ следует вычислять путем сведения их к N/4-точечным преобразованиям:
,
(39)
(40)
In this case, due to the periodicity of the sequence W nk N/4 in k with period N/4, sums (40) must be calculated only for values of k from 0 to N/4-1. Therefore, the calculation of the sequence X[k] by formulas (37), (39) and (40) requires, as it is easy to calculate, already 2N+N 2 /4 operations of multiplication and addition.
Following this path, the amount of calculations X[k] can be more and more reduced. After m=log 2 N expansions, we arrive at two-point Fourier transforms of the form
(41)
where the "single-point transformations" X 1 are simply signal samples x[n]:
X 1 = x[q]/N, q=0,1,...,N-1. (42)
As a result, we can write the FFT algorithm, which for obvious reasons is called time thinning algorithm :
X 2 \u003d (x[p] + W k 2 x) / N,
where k=0.1, p=0.1,...,N/2 -1;
X 2N/M =X N/M + W k 2N/M X N/M ,
where k=0.1,...,2N/M -1, p=0.1,...,M/2 -1;
X[k] = X N [k] =X N/2 + W k N X N/2 , (43)
where k=0,1,...,N-1
At each stage of calculations, N complex multiplications and additions are performed. And since the number of decompositions of the original sequence into half-length subsequences is log 2 N, then the total number of multiplication-add operations in the FFT algorithm is Nlog 2 N. For large N, there is a significant saving in computational operations compared to direct calculation of the DFT. For example, at N = 2 10 = 1024 the number of operations decreases by 117 times.
The FFT algorithm with time decimation considered by us is based on the calculation of the Fourier transform by forming subsequences of the input sequence x[n]. However, one can also use the subsequence decomposition of the Fourier transform X[k]. The FFT algorithm based on this procedure is called the FFT algorithm. decimation in frequency. You can read more about the fast Fourier transform, for example, in.
Random Processes and Power Spectral Density
Discrete random process x can be considered as some set or ensemble of real or complex discrete temporal (or spatial) sequences, each of which could be observed as a result of some experiment (n is the time index, i is the observation number). The sequence obtained as a result of one of the observations will be denoted by x[n]. The ensemble averaging operation (i.e. statistical averaging) will be denoted by the operator<>. Thus,
A random process is called stationary in broad sense, if its average value is constant (does not depend on time), and autocorrelation depends only on the difference of time indices m=n1-n2 (time shift or delay between samples). Thus, a broadly stationary discrete random process x[n] is characterized by a constant mean value
r xx [m] =< xx*[n] >. (44)
Note the following properties of the ACP:
r xx |r xx [m]| , r xx [-m] = r* xx [m] , (45)
which are valid for all m.
The power spectral density (PSD) is defined as the discrete-time Fourier transform (DTFT) of the autocorrelation sequence
.
(46)
The PSD, whose width is assumed to be limited to ±1/2T Hz, is a periodic function of frequency with a period of 1/T Hz. The PSD function describes the frequency distribution of the power of a random process. To confirm the name chosen for it, consider the inverse DTFT
(47)
calculated at m=0
(48)
Autocorrelation at zero shift characterizes average power random process. According to (48), the area under the curve P xx (f) characterizes the average power, therefore P xx (f) is a function of density (power per unit of frequency measurement), which characterizes the distribution of power over frequency. The pair of transformations (46) and (47) is often called Wiener-Khinchin theorem for the case of discrete time. Since r xx [-m]=r* xx [m], the PSD must be a strictly real positive function. If the AFC is a strictly real function, then r xx [-m]=r xx [m] and the PSD can be written in the form of the Fourier cosine transform
,
which also means that P xx (f) = P xx (-f), i.e. SPM is an even function.
Until now, in determining the mean value, correlation, and power spectral density of a random process, we have used statistical averaging over the ensemble. However, in practice, it is usually not possible to obtain an ensemble of implementations of the required process, from which these statistical characteristics could be calculated. It is desirable to estimate all statistical properties from one sample realization x(t), replacing y ensemble average by time average. The property that allows such a replacement to be carried out is called ergodicity. A random process is said to be ergodic if, with a probability equal to one, all its statistical characteristics can be predicted from one implementation from the ensemble using time averaging. In other words, the average values over time of almost all possible realizations of the process converge with a probability of one to the same constant value - the average value over the ensemble
. (49)
This limit, if it exists, converges to the true mean if and only if the variance of the time mean goes to zero, which means that the following condition is met:
. (50)
Here c xx [m] is the true value of the covariance of the process x[n].
Similarly, observing the value of the product of x[n] process samples at two points in time, we can expect that the average value will be equal to
(51)
The assumption of ergodicity allows not only to introduce, through time averaging, definitions for the mean and autocorrelation, but also to give a similar definition for the power spectral density
. (52)
This equivalent form of PSD is obtained by statistically averaging the modulo DTFT of the weighted data set divided by the length of the data record for the case where the number of samples increases to infinity. Statistical averaging is necessary here because the DTFT itself is a random variable that changes for each implementation of x[n]. In order to show that (52) is equivalent to the Wiener-Khinchin theorem, we represent the square of the DTFT modulus as a product of two series and change the order of the summation and statistical averaging operations:
(53)
Using the well-known expression
, (54)
relation (53) can be reduced to the following:
(55)
Note that at the last stage of the derivation of (55) we used the assumption that the autocorrelation sequence 'decays', so that
.
(56)
The relationship between the two definitions of SPM (46) and (52) is clearly shown in the diagram shown in Figure 4.
If expression (52) does not take into account the operation of mathematical expectation, then we obtain the estimate of the PSD
, (57)
which is called selective spectrum.
Rice. 4. Relationship between two methods for estimating the power spectral density
Periodogram method of spectral estimation
Above, we introduced two formal equivalent methods for determining the power spectral density (PSD). The indirect method is based on the use of an infinite data sequence to calculate an autocorrelation sequence whose Fourier transform gives the desired PSD. The direct method for determining the PSD is based on calculating the square of the modulus of the Fourier transform for an infinite sequence of data using appropriate statistical averaging. The PSD obtained without such averaging turns out to be unsatisfactory, since the root-mean-square error of such an estimate is comparable to its average value. Now we will consider averaging methods that provide smooth and statistically stable spectral estimates over a finite number of samples. JMP estimates based on direct conversion data and subsequent averaging, are called periodograms. JMP estimates, for which correlation estimates are first formed from the initial data, are called correlogram. When using any PSD estimation method, the user has to make many trade-off decisions in order to obtain statistically stable spectral estimates with the highest possible resolution from a finite number of samples. These trade-offs include, in particular, the selection of a window for weighting the data and correlation estimates and such averaging parameters in the time and frequency domains that balance the requirements for side-lobe reduction due to weighting, performing effective averaging, and providing acceptable spectral resolution. On fig. 5 is a diagram showing the main stages periodogram method
Rice. 5. The main stages of PSD estimation using the periodogram method
The application of the method begins with the collection of N data samples, which are taken at intervals of T seconds per sample, followed by (optionally) a detrending step. In order to obtain a statistically stable spectral estimate, the available data must be divided into overlapping (if possible) segments and, subsequently, the sample spectra obtained for each such segment should be averaged. The parameters of this averaging are changed by appropriate selection of the number of samples per segment (NSAMP) and the number of samples to shift the beginning of the next segment (NSHIFT), see fig. 6. The number of segments is selected depending on the required degree of smoothness (dispersion) of the spectral estimate and the required spectral resolution. With a small value of the NSAMP parameter, there are more segments to be averaged over, and therefore estimates with less variance, but also less frequency resolution, will be obtained. Increasing the segment length (NSAMP parameter) increases the resolution, naturally at the expense of increasing the estimation variance due to fewer averages. The back arrow in Fig. 5 indicates the need for several repeated passes through the data at different lengths and numbers of segments, which allows you to get more information about the process under study.
Fig.6. Splitting Data into Segments to Calculate a Periodogram
Window
One of the important issues that is common to all classical methods of spectral estimation is related to data weighting. Windowing is used to control the effects of sidelobes in spectral estimates. Note that it is convenient to consider the available finite data record as some part of the corresponding infinite sequence, visible through the applied window. So the sequence of observed data x 0 [n] from N samples can be mathematically written as the product of an infinite sequence x[n] and a rectangular window function
X 0 [n]=x[n]rect[n].
This assumes the obvious assumption that all unobserved counts are zero, regardless of whether this is actually the case. The discrete-time Fourier transform of a weighted sequence is equal to the convolution of the x[n] sequence transformations and the rectangular window rect[n]
X 0 (f)=X(f)*D N (f) , where
D N (f)= Texp(-j2pfT)sin(pfTN)/sin(pfT).
The function D N (f), called the discrete sinc function, or the Dirichlet kernel, is the DTFT of a rectangular function. The observable finite sequence transformation is a distorted version of the infinite sequence transformation. The influence of a rectangular window on a discrete-time sinusoid with a frequency f 0 is illustrated in Fig.7.
Fig.7. Illustration of the displacement of the discrete-time Fourier transform due to leakage due to data weighting.: a, c - original and weighted sequences; b, d - their Fourier transforms.
It can be seen from the figure that the sharp spectral peaks of the DTFT of the infinite sinusoidal sequence have been expanded due to convolution with the window transformation. Thus, the minimum width of the spectral peaks of a window-weighted sequence is determined by the width of the main lobe of the transformation of this window and does not depend on the data. The sidelobes of the window transformation will change the amplitudes of adjacent spectral peaks (sometimes referred to as leakage). Since the DTFT is a periodic function, the overlap of side lobes from neighboring periods can lead to additional bias. Increasing the sampling frequency reduces the effect of side-lobe superposition. Similar distortions will, of course, be observed in the case of non-sinusoidal signals. Leakage leads not only to the appearance of amplitude errors in the spectra of discrete signals, but can also mask the presence of weak signals. A number of other window functions can be proposed that can reduce the level of sidelobes compared to a rectangular window. Reducing the level of the side lobes will reduce the bias of the spectral estimate, but this comes at the cost of expanding the main lobe of the window spectrum, which naturally leads to a deterioration in resolution. Therefore, here, too, some compromise must be made between the width of the main lobe and the level of the side lobes. Several parameters are used to evaluate the quality of windows. The traditional measure is the bandwidth of the main lobe at half power. The second metric is the equivalent bandwidth entered above. Two indicators are also used to evaluate the characteristics of the side lobes. The first is their maximum level, the second is the decay rate, which characterizes the speed of the side lobes decreasing as they move away from the main lobe. Table 3 provides definitions for some commonly used discrete-time window functions, and Table 4 describes their characteristics.
Table 3 Definitions of Typical N-Point Discrete Time Window Max. side lobes level, dB -31.5
. (46)
Correlogram method estimation of the PSD is simply a substitution into expression (46) of a finite sequence of values of the autocorrelation estimate ( correlograms) instead of an infinite sequence of unknown true autocorrelation values. More information about the correlogram method of spectral estimation can be found in.
Literature
1. Rabiner L., Gould B. Theory and application of digital signal processing. M.: Mir, 1978.
2. Marple Jr. S.L. Digital spectral analysis and its applications: Per. from English. -M.: Mir, 1990.
3. Goldberg L.M., Matyushkin B.D., Polyak M.N., Digital signal processing. - M.: Radio and communication, 1990.
4. Otnes R., Enokson L. Applied analysis of time series.- M.: Mir, 1982.
Video surveillance cameras are widely used to control the traffic situation on highways with high traffic intensity. The information coming from the video cameras contains data on the temporary change in the spatial position of vehicles in the field of view of the system. The processing of this information on the basis of algorithms used in television measuring systems (TIS) makes it possible to determine the speed of vehicles and ensure traffic flow control. It is these factors that explain the growing interest in television monitoring of transport routes.
To develop methods for filtering images of vehicles against the background of noise, it is necessary to know their main parameters and characteristics. Previously, the authors carried out a study of the Fourier and wavelet spectra of natural and urban backgrounds. This work is devoted to the study of similar spectra of vehicles.
- a bank of original .bmp files of monochrome images of vehicles was created using a digital camera various types(cars and trucks, buses, for each group the number of images was 20-40 at different angles and lighting conditions); the images were 400 pixels horizontally and 300 pixels vertically; brightness range from 0 to 255 units;
- since the images contained, in addition to the vehicle, also a background component, in order to prevent its influence on the result, it was artificially suppressed to a zero level;
- the characteristics of vehicle images were analyzed by Fourier and wavelet analysis methods.
Developed in the environment MATLAB program allows you to calculate the average brightness (i.e. the mathematical expectation of image brightness), brightness dispersion, Fourier spectrum of individual and total image lines, spectrograms, as well as wavelet spectra using various well-known wavelets (Haar, Daubechies, Simlet, etc.). The results of the analysis are reflected in the form of two-dimensional and 3D image spectra.
Based on the research results, the following conclusions can be drawn:
- average brightness characteristics (average brightness, dispersion) of images of different vehicles have similar values for all types; a significant impact on the brightness characteristics have solar glare from the windows and surfaces of the car; depending on the intensity and direction of lighting, black cars may have brightness characteristics similar to light cars;
- regardless of the type of vehicle, Fourier and wavelet spectra have a similar structure;
- the Fourier width of the vehicle spectrum weakly depends on the vehicle type; the spectrum has a significantly non-uniform structure that changes with changes in lighting and vehicle orientation; the spectrum in the horizontal plane has a more uneven structure than in the vertical one; the spectral characteristics of semi-trucks and buses are greatly influenced by drawings and inscriptions (advertising) on its surfaces;
- when turning cars, the change in the spectra of images in the horizontal plane is significant, the spectrum in the vertical plane remains quite stable; this is especially well seen in the wavelet spectra;
- analysis of the spectra of an individual vehicle and a vehicle against the background of interference shows that they differ in the amplitude levels of the spectral components; in the absence of a background, the vertical spectrum is much more uniform; for images of cars without a background, there is a greater probability of deep dips in the spectrum (higher unevenness), the envelope of the spectrum of images with a background is more uniform than without a background;
- studies have shown that due to strong influence a large number of factors, the spectral characteristics of vehicles (both obtained using Fourier analysis and wavelet analysis) do not allow us to identify stable spectral features of vehicle images; this reduces the efficiency of the spectral filtering of images, carried out to suppress the background;
- V automated systems For traffic control, in order to distinguish cars against the background of interference, it is necessary to use a set of features, such as color, spectrum, geometric parameters of objects (sizes and ratios of sizes) and dynamic characteristics.
BIBLIOGRAPHY
- Makaretsky E.A., Nguyen L.Kh. Study of the characteristics of images of natural and urban backgrounds / / Izv. Tulsk. State. University. Radio engineering and radio optics. - Tula, 2005. - T. 7.- P. 97-104.
Bibliographic link
Makaretsky E.A. STUDY OF THE FOURIER AND WAVELET SPECTRA OF IMAGES OF VEHICLES // Basic Research. - 2006. - No. 12. - P. 80-81;URL: http://fundamental-research.ru/ru/article/view?id=5557 (date of access: 01/15/2020). We bring to your attention the journals published by the publishing house "Academy of Natural History"
Spectral analysis
Spectral analysis is a wide class of data processing methods based on their frequency representation, or spectrum. The spectrum is obtained by decomposing the original function, depending on time (time series) or spatial coordinates (for example, images), into the basis of some periodic function. The most commonly used for spectral processing is the Fourier spectrum obtained on the basis of the sine basis (Fourier expansion, Fourier transform).
The main meaning of the Fourier transform is that the original non-periodic function of an arbitrary form, which cannot be described analytically and therefore is difficult to process and analyze, is represented as a set of sines or cosines with different frequencies, amplitudes and initial phases.
In other words, a complex function is transformed into a set of simpler ones. Each sinusoid (or cosine wave) with a certain frequency and amplitude, obtained as a result of the Fourier decomposition, is called spectral component or harmonica. The spectral components form Fourier spectrum.
Visually, the Fourier spectrum is represented as a graph, on which the circular frequency, denoted by the Greek letter "omega", is plotted along the horizontal axis, and the amplitude of the spectral components, usually denoted by the Latin letter A, is plotted along the vertical axis. Then each spectral component can be represented as a reference, position which horizontally corresponds to its frequency, and its height corresponds to its amplitude. A harmonic with zero frequency is called constant component(in time representation it is a straight line).
Even a simple visual analysis of the spectrum can tell a lot about the nature of the function from which it was derived. It is intuitively clear that rapid changes in the initial data give rise to components in the spectrum with high frequency, and slow ones with low. Therefore, if the amplitude of the components in it rapidly decreases with increasing frequency, then the original function (for example, a time series) is smooth, and if the spectrum contains high-frequency components with a large amplitude, then the original function will contain sharp fluctuations. So, for a time series, this may indicate a large random component, the instability of the processes described by it, the presence of noise in the data.
Spectral manipulation is based on spectrum manipulation. Indeed, if we reduce (suppress) the amplitude of the high-frequency components, and then, based on the modified spectrum, restore the original function by performing the inverse Fourier transform, then it will become smoother due to the removal of the high-frequency component.
For a time series, for example, this means removing information about daily sales, which are highly affected by random factors, and leaving more stable trends, such as seasonality. You can, on the contrary, suppress components with a low frequency, which will allow you to remove slow changes, and leave only fast ones. In the case of a time series, this would mean the suppression of the seasonal component.
By applying the spectrum in this way, the desired change in the original data can be achieved. The most commonly used time series smoothing is by removing or reducing the amplitude of high-frequency components in the spectrum.
To manipulate spectra, filters are used - algorithms that can control the shape of the spectrum, suppress or enhance its components. chief property any filter is its amplitude-frequency characteristic (AFC), the shape of which depends on the transformation of the spectrum.
If the filter passes only spectral components with a frequency below a certain cutoff frequency, then it is called a low-pass filter (LPF), and it can be used to smooth data, clean it from noise and anomalous values.
If the filter passes the spectral components above a certain cutoff frequency, then it is called a high-pass filter (HPF). It can be used to suppress slow changes, such as seasonality in a data series.
In addition, many other types of filters are used: mid-pass filters, trap filters and band-pass filters, as well as more complex ones that are used in signal processing in electronics. By selecting the type and shape of the filter's frequency response, it is possible to achieve the desired transformation of the original data through spectral processing.
When performing frequency filtering of data in order to smooth and remove noise, it is necessary to correctly specify the bandwidth of the low-pass filter. If it is set too high, then the degree of smoothing will be insufficient, and the noise will not be completely suppressed. If it is too narrow, then along with the noise, changes that bring useful information. If in technical applications There are strict criteria for determining the optimal characteristics of filters, then in analytical technologies it is necessary to use mainly experimental methods.
Spectral analysis is one of the most efficient and well developed methods of data processing. Frequency filtering is just one of its many applications. In addition, it is used in correlation and statistical analysis, signal and function synthesis, model building, etc.
1. Fourier transform and signal spectrum
In many cases, the task of obtaining (calculating) the signal spectrum is as follows. There is an ADC, which with a sampling frequency Fd converts a continuous signal arriving at its input during the time T, into digital readings - N pieces. Next, the array of readings is fed into a certain program that gives out N / 2 of some numerical values (the programmer who pulled from internet wrote a program, claims that it does the Fourier transform).To check if the program is working correctly, we will form an array of readings as the sum of two sinusoids sin(10*2*pi*x)+0.5*sin(5*2*pi*x) and slip it into the program. The program drew the following:
fig.1 Graph of the time function of the signal
fig.2 Graph of signal spectrum
On the spectrum graph there are two sticks (harmonics) 5 Hz with an amplitude of 0.5 V and 10 Hz - with an amplitude of 1 V, all as in the formula of the original signal. Everything is fine, well done programmer! The program is working correctly.
This means that if we apply a real signal from a mixture of two sinusoids to the input of the ADC, then we will get a similar spectrum consisting of two harmonics.
Total, our real measured signal, duration 5 sec, digitized by the ADC, i.e. represented discrete counts, has discrete non-periodic range.
From a mathematical point of view, how many mistakes are there in this phrase?
Now the authorities have decided we decided that 5 seconds is too long, let's measure the signal in 0.5 seconds.
fig.3 Graph of the function sin(10*2*pi*x)+0.5*sin(5*2*pi*x) for a measurement period of 0.5 sec
fig.4 Function spectrum
Something is not right! The 10 Hz harmonic is drawn normally, but instead of a 5 Hz stick, several incomprehensible harmonics appeared. We look on the Internet, what and how ...
In, they say that zeros must be added to the end of the sample and the spectrum will be drawn normal.
fig.5 Finished zeros up to 5 seconds
fig.6 We got the spectrum
Still not what it was at 5 seconds. You have to deal with the theory. Let's go to Wikipedia- source of knowledge.
2. A continuous function and its representation by a Fourier series
Mathematically, our signal with a duration of T seconds is a certain function f(x) given on the interval (0, T) (X in this case is time). Such a function can always be represented as a sum of harmonic functions (sine or cosine) of the form:
(1), where:
K - number of trigonometric function (number of harmonic component, harmonic number)
T - segment where the function is defined (signal duration)
Ak - amplitude of k-th harmonic component,
?k - initial phase of k-th harmonic component
What does it mean to "represent a function as a sum of a series"? This means that by adding the values of the harmonic components of the Fourier series at each point, we will get the value of our function at this point.
(More strictly, the standard deviation of the series from the function f(x) will tend to zero, but despite the standard convergence, the Fourier series of the function, generally speaking, is not required to converge pointwise to it. See https://ru.wikipedia.org/ wiki/Fourier_Series .)
This series can also be written as:
(2),
Where , k-th complex amplitude.
The relationship between the coefficients (1) and (3) is expressed by the following formulas:
Note that all these three representations of the Fourier series are completely equivalent. Sometimes, when working with Fourier series, it is more convenient to use the exponents of the imaginary argument instead of sines and cosines, that is, to use the Fourier transform in complex form. But it is convenient for us to use formula (1), where the Fourier series is represented as a sum of cosine waves with the corresponding amplitudes and phases. In any case, it is incorrect to say that the result of the Fourier transform of the real signal will be the complex amplitudes of the harmonics. As the wiki correctly says, "The Fourier transform (?) is an operation that maps one function of a real variable to another function, also of a real variable."
Total:
The mathematical basis of the spectral analysis of signals is the Fourier transform.
The Fourier transform allows us to represent a continuous function f(x) (signal) defined on the segment (0, T) as the sum of an infinite number (infinite series) of trigonometric functions (sine and/or cosine) with certain amplitudes and phases, also considered on the segment (0, T). Such a series is called a Fourier series.
We note some more points, the understanding of which is required for the correct application of the Fourier transform to signal analysis. If we consider the Fourier series (the sum of sinusoids) on the entire X-axis, then we can see that outside the segment (0, T), the function represented by the Fourier series will periodically repeat our function.
For example, in the graph in Fig. 7, the original function is defined on the segment (-T \ 2, + T \ 2), and the Fourier series represents a periodic function defined on the entire x-axis.
This is because the sinusoids themselves are periodic functions, respectively, and their sum will be a periodic function.
fig.7 Representation of a non-periodic original function by a Fourier series
Thus:
Our original function is continuous, non-periodic, defined on some interval of length T.
The spectrum of this function is discrete, that is, it is presented as an infinite series of harmonic components - the Fourier series.
In fact, a certain periodic function is defined by the Fourier series, which coincides with ours on the segment (0, T), but this periodicity is not essential for us.
The periods of the harmonic components are multiples of the segment (0, T) on which the original function f(x) is defined. In other words, the harmonic periods are multiples of the duration of the signal measurement. For example, the period of the first harmonic of the Fourier series is equal to the interval T on which the function f(x) is defined. The period of the second harmonic of the Fourier series is equal to the interval T/2. And so on (see Fig. 8).
fig.8 Periods (frequencies) of the harmonic components of the Fourier series (here T = 2?)
Accordingly, the frequencies of the harmonic components are multiples of 1/T. That is, the frequencies of the harmonic components Fk are equal to Fk= k\T, where k ranges from 0 to?, for example, k=0 F0=0; k=1 F1=1\T; k=2 F2=2\T; k=3 F3=3\T;… Fk= k\T (at zero frequency - constant component).
Let our original function be a signal recorded for T=1 sec. Then the period of the first harmonic will be equal to the duration of our signal T1=T=1 sec and the frequency of the harmonic is 1 Hz. The period of the second harmonic will be equal to the duration of the signal divided by 2 (T2=T/2=0.5 sec) and the frequency is 2 Hz. For the third harmonic T3=T/3 sec and the frequency is 3 Hz. And so on.
The step between harmonics in this case is 1 Hz.
Thus, a signal with a duration of 1 sec can be decomposed into harmonic components (to obtain a spectrum) with a frequency resolution of 1 Hz.
To increase the resolution by 2 times to 0.5 Hz, it is necessary to increase the measurement duration by 2 times - up to 2 seconds. A signal with a duration of 10 seconds can be decomposed into harmonic components (to obtain a spectrum) with a frequency resolution of 0.1 Hz. There are no other ways to increase the frequency resolution.
There is a way to artificially increase the duration of the signal by adding zeros to the array of samples. But it does not increase the real frequency resolution.
3. Discrete signals and discrete Fourier transform
With the development of digital technology, the ways of storing measurement data (signals) have also changed. If earlier the signal could be recorded on a tape recorder and stored on tape in analog form, now the signals are digitized and stored in files in the computer's memory as a set of numbers (counts).The usual scheme for measuring and digitizing a signal is as follows.
fig.9 Scheme of the measuring channel
The signal from the measuring transducer arrives at the ADC during a period of time T. The signal samples (sample) obtained during the time T are transferred to the computer and stored in memory.
fig.10 Digitized signal - N readings received in time T
What are the requirements for signal digitization parameters? A device that converts an input analog signal into a discrete code (digital signal) is called an analog-to-digital converter (ADC, English Analog-to-digital converter, ADC) (Wiki).
One of the main parameters of the ADC is the maximum sampling rate (or sampling rate, English sample rate) - the frequency of taking samples of a signal continuous in time during its sampling. Measured in hertz. ((Wiki))
According to the Kotelnikov theorem, if a continuous signal has a spectrum limited by the frequency Fmax, then it can be completely and uniquely restored from its discrete samples taken at time intervals 
, i.e. with frequency Fd ? 2*Fmax, where Fd - sampling frequency; Fmax - maximum frequency of the signal spectrum. In other words, the signal sampling rate (ADC sampling rate) must be at least 2 times the maximum frequency of the signal that we want to measure.
And what will happen if we take readings with a lower frequency than required by the Kotelnikov theorem?
In this case, the effect of "aliasing" (aka stroboscopic effect, moiré effect) occurs, in which the high-frequency signal after digitization turns into a low-frequency signal that does not actually exist. On fig. 5 high frequency red sine wave is the real signal. The lower frequency blue sine wave is a dummy signal resulting from the fact that more than half a period of a high-frequency signal has time to pass during the sampling time.
Rice. 11. The appearance of a false low-frequency signal when the sampling rate is not high enough
To avoid the effect of aliasing, a special anti-aliasing filter is placed in front of the ADC - LPF (low-pass filter), which passes frequencies below half the ADC sampling frequency, and cuts off higher frequencies.
In order to calculate the spectrum of a signal from its discrete samples, the discrete Fourier transform (DFT) is used. We note once again that the spectrum of a discrete signal is "by definition" limited by the frequency Fmax, which is less than half the sampling frequency Fd. Therefore, the spectrum of a discrete signal can be represented by the sum of a finite number of harmonics, in contrast to the infinite sum for the Fourier series of a continuous signal, the spectrum of which can be unlimited. According to the Kotelnikov theorem, the maximum harmonic frequency must be such that it accounts for at least two samples, so the number of harmonics is equal to half the number of samples of the discrete signal. That is, if there are N samples in the sample, then the number of harmonics in the spectrum will be equal to N/2.
Consider now the discrete Fourier transform (DFT).
Comparing with the Fourier series
We see that they coincide, except that the time in the DFT is discrete and the number of harmonics is limited to N/2 - half the number of samples.
The DFT formulas are written in dimensionless integer variables k, s, where k are the numbers of signal samples, s are the numbers of spectral components.
The value of s shows the number of full oscillations of the harmonic in the period T (the duration of the signal measurement). The discrete Fourier transform is used to find the amplitudes and phases of harmonics numerically, i.e. "on the computer"
Returning to the results obtained at the beginning. As mentioned above, when expanding a non-periodic function (our signal) into a Fourier series, the resulting Fourier series actually corresponds to a periodic function with period T. (Fig. 12).
fig.12 Periodic function f(x) with period Т0, with measurement period Т>T0
As can be seen in Fig. 12, the function f(x) is periodic with period Т0. However, due to the fact that the duration of the measurement sample T does not coincide with the period of the function T0, the function obtained as a Fourier series has a discontinuity at the point T. As a result, the spectrum of this function will contain a large number of high frequency harmonics. If the duration of the measurement sample T coincided with the period of the function T0, then only the first harmonic (a sinusoid with a period equal to the sample duration) would be present in the spectrum obtained after the Fourier transform, since the function f(x) is a sinusoid.
In other words, the DFT program "does not know" that our signal is a "piece of a sine wave", but is trying to represent a periodic function as a series, which has a gap due to the inconsistency of the individual pieces of the sine wave.
As a result, harmonics appear in the spectrum, which in total should represent the form of the function, including this discontinuity.
Thus, in order to obtain the "correct" spectrum of the signal, which is the sum of several sinusoids with different periods, it is necessary that an integer number of periods of each sinusoid fit on the signal measurement period. In practice, this condition can be met for a sufficiently long duration of the signal measurement.
Fig.13 An example of the function and spectrum of the signal of the kinematic error of the gearbox
With a shorter duration, the picture will look "worse":
Fig.14 An example of the function and spectrum of the rotor vibration signal
In practice, it can be difficult to understand where are the “real components” and where are the “artifacts” caused by the non-multiplicity of the periods of the components and the duration of the signal sample or the “jumps and breaks” of the waveform. Of course, the words "real components" and "artifacts" are not in vain quoted. The presence of many harmonics on the spectrum graph does not mean that our signal actually “consists” of them. It's like thinking that the number 7 "consists" of the numbers 3 and 4. The number 7 can be represented as the sum of the numbers 3 and 4 - this is correct.
So is our signal ... or rather, not even “our signal”, but a periodic function compiled by repeating our signal (sampling) can be represented as a sum of harmonics (sinusoids) with certain amplitudes and phases. But in many cases important for practice (see the figures above), it is indeed possible to relate the harmonics obtained in the spectrum to real processes that are cyclic in nature and make a significant contribution to the signal shape.
Some results
1. The real measured signal, duration T sec, digitized by the ADC, that is, represented by a set of discrete samples (N pieces), has a discrete non-periodic spectrum, represented by a set of harmonics (N/2 pieces).2. The signal is represented by a set of real values and its spectrum is represented by a set of real values. The harmonic frequencies are positive. The fact that it is more convenient for mathematicians to represent the spectrum in a complex form using negative frequencies does not mean that “it is right” and “it should always be done this way”.
3. The signal measured on the time interval T is determined only on the time interval T. What happened before we began to measure the signal, and what will happen after that - this is unknown to science. And in our case - it is not interesting. The DFT of a time-limited signal gives its "real" spectrum, in the sense that, under certain conditions, it allows you to calculate the amplitude and frequency of its components.
Used materials and other useful materials.
