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 heat 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 described 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, a change in the orbit of a comet is due to the influence of the attraction of several planets, a change in the flow of heat is 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 summed up to give a solution to the 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 of complex Fourier methods, scientists have been 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. In the future, these methods were used 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.
The Introductory Overview section discusses two very simple examples (taken from Shumway, 1988) to illustrate the nature of spectral analysis and the interpretation of the results. If you are unfamiliar with this method, it is recommended that you review this section of this chapter first.
Review and data file. The Sunspot.sta file contains a fraction of the known sunspot numbers (Wolfer) from 1749 to 1924 (Anderson, 1971). Below is a list of the first few data from the example file.
It is assumed that the number of sunspots affects the weather on earth, as well as agriculture, telecommunications, etc. Using this analysis, one can try to find out whether sunspot activity is indeed cyclic in nature (in fact it is, these data are widely discussed in the literature; see, for example, Bloomfield, 1976, or Shumway, 1988).
Analysis definition. After running the analysis, open the Sunspot.sta data file. Click the Variables button and select the Spots variable (note that if the data file Sunspot.sta is the currently open data file, and the Spots variable is the only variable in that file, Spots will automatically be selected when the Time Series Analysis dialog box opens). Now click on the Fourier (spectral) analysis button to open the Fourier (spectral) analysis dialog box.
Before applying spectral analysis, first plot the number of sunspots. Note that the Sunspot.sta file contains the corresponding years as observation names. To use these names in line plots, click the View Series tab and select Case Names under Label Points. Also, select Set x-axis scale manually and Min. = 1 and Step = 10. Then click the Graph button next to the Preview highlight button. variable.
The number of sunspots seems to follow a cyclic pattern. There is no trend, so go back to the Spectrum Analysis window and deselect the Remove Linear Trend option in the Transform Original Series group.
Obviously, the mean of the series is greater than 0 (zero). Therefore, leave the Subtract mean option selected [otherwise the periodogram will "clog" with a very large peak at frequency 0 (zero)].
Now you are ready to start the analysis. Now click OK (One Dimensional Fourier Analysis) to bring up the Fourier Spectral Analysis Results dialog box.
View results. The information section at the top of the dialog box shows some summary statistics for the series. It also shows the five largest periodogram peaks (by frequency). The largest three peaks are at frequencies 0.0852, 0.0909 and 0.0114. This information is often useful when analyzing very large series (for example, those with more than 100,000 observations) that are not easily plotted in one plot. In this case, however, it is easy to see the periodogram values; by clicking the Periodogram button under Periodogram and Spectral Density Plots.
The periodogram plot shows two distinct peaks. The maximum is at a frequency of approximately 0.9. Return to the Spectral Analysis Results window and click the Summary button to see all periodogram values (and other results) in the results table. Shown below is a portion of the results table with the largest peak set from the periodogram.
As discussed in the Introductory Overview section, Frequency is the number of cycles per unit of time (where each observation is one unit of time). Thus, a Frequency of 0.0909 corresponds to a value of 11 Period (the number of units of time required for a complete cycle). Since the sunspot data in Sunspot.sta are annual observations, it can be concluded that there is a pronounced 11-year (perhaps slightly longer than 11-year) cycle in sunspot activity.
Spectral density. Usually, to calculate spectral density estimates, the periodogram is smoothed to remove random fluctuations. The weighted moving average type and window width can be selected in the Spectral Windows section. The Introductory Overview section discusses these options in detail. For our example, let's leave the default window selected (Hamming width 5) and select the Spectral Density plot.
The two peaks are now even clearer. Let's look at the values of the periodogram over the period. Highlight the Period field in the Graph section. Now select the Spectral Density plot.
Again, there is a pronounced 11-year cycle in sunspot activity; moreover, there are signs of a longer cycle of about 80-90 years.
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 law of conservation of energy 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 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 low-pass filter) 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 region (for example, time) lead to a periodic continuation in the transformation region (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 the time and frequency domains, we obtain
. (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 NTPF 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 in fact 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 doubling 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 of an 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 an TDGF 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 equal to 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 - time index, i - 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 change to be made 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 DTFT modulus 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. PSD estimates based on direct data transformation 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 compromise 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 efficient 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.
Fourier transform- this is a family of mathematical methods based on the decomposition of the original continuous function of time into a set of basic harmonic functions (which are sinusoidal functions) of various frequencies, amplitudes and phases. It can be seen from the definition that the main idea of the transformation is that any function can be represented as an infinite sum of sinusoids, each of which will be characterized by its amplitude, frequency and initial phase.
The Fourier transform is the founder of spectral analysis. Spectral analysis is a signal processing method that allows you to characterize the frequency content of the measured signal. Depending on how the signal is represented, different Fourier transforms are used. There are several types of Fourier transform:
– Continuous Fourier Transform (in English literature Continue Time Fourier Transform – CTFT or, for short, FT);
– Discrete Fourier Transform (in English literature Discrete Fourier Transform – DFT);
– Fast Fourier transform (in English literature Fast Fourier transform – FFT).
Continuous Fourier Transform
The Fourier transform is a mathematical tool used in various scientific fields. In some cases, it can be used as a means of solving complex equations that describe dynamic processes that occur under the influence of electrical, thermal or light energy. In other cases, it allows you to highlight the regular components in a complex oscillatory signal, so that you can correctly interpret experimental observations in astronomy, medicine and chemistry. A continuous transformation is actually a generalization of Fourier series, provided that the period of the expanded function tends to infinity. Thus, the classical Fourier transform deals with the spectrum of the signal taken over the entire range of the existence of the variable.
There are several types of writing the continuous Fourier transform, which differ from each other by the value of the coefficient in front of the integral (two forms of writing):
or 
where and is the Fourier image of the function or the frequency spectrum of the function ;
- circular frequency.
It should be noted that different types of recording are found in various fields of science and technology. The normalization factor is necessary for the correct scaling of the signal from the frequency domain to the time domain. The normalization factor reduces the amplitude of the signal at the output of the inverse transformation so that it matches the amplitude of the original signal. In the mathematical literature, the direct and inverse Fourier transforms are multiplied by the factor , while in physics, most often, the factor is not set for the direct transformation, but the factor is set for the reverse. If we sequentially calculate the direct Fourier transform of a certain signal, and then take the inverse Fourier transform, then the result of the inverse transform should completely coincide with the original signal.
If the function is odd on the interval (−∞, +∞), then the Fourier transform can be represented in terms of the sine function:
If the function is even on the interval (−∞, +∞), then the Fourier transform can be represented in terms of the cosine function:
Thus, the continuous Fourier transform allows us to represent a non-periodic function as an integral of a function representing at each of its points the coefficient of the Fourier series for a non-periodic function.
The Fourier transform is reversible, that is, if its Fourier image was calculated from the function, then the original function can be uniquely restored from the Fourier image. The inverse Fourier transform is understood as an integral of the form (two forms of writing):
or 
where is the Fourier image of the function or the frequency spectrum of the function ;
- circular frequency.
If the function is odd on the interval (−∞, +∞), then the inverse Fourier transform can be represented in terms of the sine function:
If the function is even on the interval (−∞, +∞), then the inverse Fourier transform can be represented in terms of the cosine function:
As an example, consider the following function 
. The graph of the exponential function under study is presented below.
Since the function is an even function, then the continuous Fourier transform will be defined as follows:
As a result, we obtained the dependence of the change in the studied exponential function on the frequency interval (see below).
The continuous Fourier transform is usually used in theory when considering signals that change in accordance with given functions, but in practice they usually deal with measurement results that represent discrete data. The measurement results are recorded at regular intervals with a certain sampling frequency, for example, 16000 Hz or 22000 Hz. However, in the general case, discrete readings can go unevenly, but this complicates the mathematical apparatus of analysis, so it is usually not used in practice.
There is an important theorem of Kotelnikov (in the foreign literature there is the name “Nyquist-Shannon theorem”, “sample theorem”), which states that an analog periodic signal with a finite (limited in width) spectrum (0 ... fmax) can be uniquely restored without distortions and losses in their discrete readings, taken with a frequency greater than or equal to twice the upper frequency of the spectrum - sampling frequency (fdisc >= 2*fmax). In other words, at a sampling rate of 1000 Hz, a signal with a frequency of up to 500 Hz can be recovered from an analog periodic signal. It should be noted that the discretization of a function in time leads to the periodization of its spectrum, and the discretization of the spectrum in frequency leads to the periodization of the function.
This is one of the Fourier transforms widely used in digital signal processing algorithms.
The direct discrete Fourier transform associates a time function , which is defined by N-measurement points on a given time interval, with another function , which is defined on a frequency interval. It should be noted that the function on the time interval is specified using N-samples, and the function on the frequency domain is specified using the K-fold spectrum.
k ˗ frequency index.
The frequency of the kth signal is determined by the expression
where T is the period of time during which the input data were taken.
The direct discrete transform can be rewritten in terms of the real and imaginary components. The real component is an array containing the values of the cosine components, and the imaginary component is an array containing the values of the sine components.
It can be seen from the last expressions that the conversion decomposes the signal into sinusoidal components (called harmonics) with frequencies from one oscillation per period to N oscillations per period.
The discrete Fourier transform has a feature, since a discrete sequence can be obtained by the sum of functions with different composition of the harmonic signal. In other words, a discrete sequence is decomposed into harmonic variables - ambiguously. Therefore, when decomposing a discrete function using a discrete Fourier transform, high-frequency components appear in the second half of the spectrum, which were not in the original signal. This high-frequency spectrum is a mirror image of the first part of the spectrum (in terms of frequency, phase and amplitude). Usually the second half of the spectrum is not considered, and the signal amplitudes of the first part of the spectrum are doubled.
It should be noted that the expansion of a continuous function does not lead to the appearance of a mirror effect, since a continuous function is uniquely decomposed into harmonic variables.
The amplitude of the DC component is the average value of the function over a selected period of time and is determined as follows:
The amplitudes and phases of the frequency components of the signal are determined by the following relationships:
The resulting amplitude and phase values are called polar notation. The resulting signal vector will be defined as follows:
Consider the algorithm for transforming a discretely given function on a given interval (on a given period) with the number of initial points
D spark Fourier transform
As a result of the transformation, we obtain the real and imaginary values of the function , which is defined on the frequency range.
The inverse discrete Fourier transform associates a frequency function , which is defined by a K-fold spectrum on the frequency domain, with another function , which is defined on the time domain.
N ˗ the number of signal values measured per period, as well as the multiplicity of the frequency spectrum;
k ˗ frequency index.
As already mentioned, the discrete Fourier transform maps N-points of a discrete signal to N-complex spectral samples of the signal. Computing one spectral sample requires N operations of complex multiplication and addition. Thus, the computational complexity of the discrete Fourier transform algorithm is quadratic, in other words, complex multiplication and addition operations are required.
1Video 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.
- using a digital camera, a bank of original .bmp files of monochrome images of vehicles of various types was created (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.
The program developed in the MATLAB environment 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, and wavelet spectra using various known wavelets (Haar, Daubechies, Simlet and 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, the probability of deep dips in the spectrum is higher (higher unevenness), the envelope of the spectrum of images with a background is more uniform than without a background;
- the conducted studies have shown that due to the strong influence of 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 spectral image filtering for background suppression;
- in automated traffic control systems, in order to distinguish vehicles against the background of interference, it is necessary to use a set of features, such as color, spectrum, geometric parameters of objects (sizes and size ratios) 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 // Fundamental 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"
