How to calculate the signal power spectrum
Using the FFT to find the real and imaginary parts of the signal spectrum, the square of the real part of the square price of the imaginary part is the power spectrum.
The spectrum of a periodic continuous signal x(t) can be expressed as a discrete non-periodic sequence Xn, its amplitude spectrum of the power spectrum │ Xn│2 square of the power spectrum of the sequence, it is called the periodic signal “power spectrum”.
The discrete spectrum of the periodic signal Xn is calculated by the Fourier transform equation,
T+t0
Xn=(1/T)∫x(t)ej(2πnt/T)dt
t0
The period of the periodic signal x(t) is indicated by T. The period of the periodic signal x(t) is called the fundamental frequency. n is the discrete frequency. “n is the independent variable of the discrete spectrum, which takes only integer values and represents a multiple of the fundamental frequency.
In general, the discrete spectrum Xn is a complex number, can be used to mode │ Xn │ and amplitude angle θn expressed as Xn = │ Xn │ ej θn, the two were called the amplitude spectrum and phase spectrum. For example, you can calculate the amplitude of 1, the average value of zero periodic square pulse amplitude spectrum for
│Xn│=2/πn,n=odd
=0, the rest of the
So its power spectrum is
│Xn│2=4/(πn)2n=odd,
=0, the rest of the
After the proposal of Fourier series, first people observe the natural world, and then they will be able to see the phase and phase spectrum of the pulse. In the late 19th century, Schuster proposed to use the magnitude square of the Fourier series as a measure of the power in a function, and named it “periodogram”. This was the earliest formulation of classical spectral estimation, and it is still used today, except that the discrete Fourier transform (DFT) is now computed using the fast Fourier transform (FFT), and the magnitude square of the DFT is used as a measure of power in the signal.
The poor variance performance of the periodogram led to the investigation of alternative methods of analysis, and in 1927 Yule proposed a linear regression equation for modeling a time series.Yule’s work actually became the basis for the most important method of modern spectral estimation, parametric modeling for spectral estimation.
Walker used Yule’s analytical method to study the decaying sinusoidal time series, resulting in the Yule-Walker equation. It can be said that both Yule and Walker were pioneers in pioneering autoregressive modeling.
Methods for computing the fast Fourier transform
Fast methods for computing the discrete Fourier transform are the FFT algorithm extracted by time and the FFT algorithm extracted by frequency. The former arranges the time-domain signal sequence by even-odd and the latter arranges the frequency-domain signal sequence by even-odd. They both draw on two characteristics: first, periodicity; second, symmetry, where the symbol * represents its conjugate. In this way, the calculation of the discrete Fourier transform can be divided into a number of steps, the efficiency of the calculation is greatly improved.
Time Extraction Algorithm Let the length of the signal sequence is N=2, where M is a positive integer, the time domain signal sequence x(n) can be decomposed into two parts, one is the even part of x(2n), and the other is the odd part of x(2n+1), and then the discrete Fourier transform of the signal sequence x(n) can be expressed and computed with the discrete Fourier transform of the two N/2 sampling points. Considering and the periodicity of the discrete Fourier transform, Eq. (1) can be written as
(3) where (4a) (4b) It follows that Eq. (4) is two discrete Fourier transforms containing only N/2 points, with G(k) including only the sequence of even points in the original signal sequence, and H(k) including only its sequence of odd points. Although k = 0, 1, 2, …, N-1, both G(k) and H(k) have a period of N/2, and their values repeat with N/2 cycles.
For then from Eqs. (3) and (4) we get (5a)(5b)
Thus, the discrete Fourier transform of a sequence of signals x(n) with N sampling points can be derived from the discrete Fourier transform of a sequence of two N/2 sampling points. By analogy, this time-dependent extraction algorithm calculates the discrete Fourier transform by dividing the input signal sequence into smaller and smaller subsequences, which are finally synthesized into the discrete Fourier transform of N points.
The discrete Fourier transform operation of Eq. (5) is usually represented by the signal flow diagram of the butterfly algorithm in Fig. 1. For example, the discrete Fourier transform of the signal sequence x(n) for a sampling point of N=8=2 can be calculated using the signal flow diagram of the FET algorithm shown in Figure 2.
① The calculation of the discrete Fourier transform of the N=2 points consists entirely of butterfly operations, and requires M levels of operations, each including N/2 butterfly operations, with a total of one butterfly operation. So, the total computation is subcomplex multiplication and Nlog2N subcomplex addition.
② The FFT algorithm proceeds iteratively by stages, and the computational formula can be written as
(6) The input signal of the N sampling points has N original data x0(n), and after the first stage of the operation, it results in a new N data x1(n), and then after the second stage of the iterative operation, it obtains another N data x2(n), and so on, until the final result of x(k) = xM(k) = X(k) in the step-by-step iterative calculation, the output data of each butterfly operation is stored in the original input data storage unit, the implementation of the so-called “i.e., bitwise calculation”, which can save a lot of intermediate data storage registers.
③ The weighting coefficients in the butterfly operation increase exponentially with the number of iteration levels. The change rule of the coefficients can be seen from Figure 2. For the N=8,M=3 case, three levels of iterative operations are required. In the first iteration, only one kind of weighting coefficient is used; the span interval of butterfly operation is equal to 1. In the second iteration, two kinds of weighting coefficients, i.e.,,; the span interval of butterfly operation is equal to 2. In the third iteration, four kinds of different weighting coefficients, i.e.,,,; the span interval of butterfly operation is equal to 4. It can be seen that the number of different weighting coefficients in each iteration is doubled compared with the previous iteration; the span interval is also doubled. interval is also doubled.
④ Input data sequence x (n) need to be rearranged as x (0), x (4), x (2), x (6), x (1), x (5), x (3), x (7), this is in accordance with the code position of the binary number of inverted inverted by the reverse number of the order, for example, N = 8 in the number of “1” in the binary number is “001 For example, the binary number of “1” in N=8 is “001”, and its code position is inverted to “100”, which is the decimal number “4”.
Frequency Extraction Algorithm The FFT algorithm of frequency extraction is to decompose the frequency domain signal sequence X(k) into odd and even parts, but the algorithm is still a step-by-step operation from the beginning of the time domain signal sequence, and the same is to split the N points into N/2 points to calculate the FFT, which can be used to directly calculate the N multiplication required by the Discrete Fourier Transform is reduced to the times.
In the case of N=2, the N-point input sequence x(n) is divided into two halves before and after
(7)
The length of the time sequence x1(n)±x2(n) is N/2, and so the discrete Fourier transform of the N-points can be written as
(8a)
(8b)
Frequency signal sequence X(2l) is time signal sequence x1(n)+x2(n) of the N/2-point discrete Fourier transform, and the frequency signal sequence X(2l+1) is the N/2-point discrete Fourier transform of the time signal sequence [x1(n)-x2(n)], so the N-point discrete Fourier transform is computed by adding (subtracting) and multiplying twice and multiplying once to obtain the two sub-sequences from the original sequence, and so the frequency extraction algorithm also has the the form of butterfly operation. The basic butterfly operation formula of FFT with base 2 is
(9)
The computation is exactly the same as that of the time extraction algorithm, i.e., it requires only times multiplication operation and Nlog2N times addition (subtraction) operation. Figure 3 represents the signal flow diagram of the discrete Fourier transform for N=8=2 points. As can be seen from the figure, it performs i.e., bitwise computation with three levels of iteration, where the input data is stored in natural order, the coefficients used are also in natural order, and the final result is stored in reverse binary order.
In fact, there is a transposition relationship between the signal flow graphs of the frequency extraction algorithm and the time extraction algorithm, and a variety of geometries can be derived if the flow graphs are appropriately deformed.
In addition to the base-2 FFT algorithm, there are also high-base FFT algorithms such as base-4, base-8, and FFT algorithms with an arbitrary number as the base.
Fast Fourier Transform FFT (FastFourierTransform)
Looked at a lot of online explanations about the FFT, some are directly ignored the derivation of the formula, and some other derivations, but the details of the derivation is not clear. With the mindset of learning without understanding, I will make the details of FFT thinking and derivation clear with formulas, so that future generations can learn FFT in more detail.
Before understanding FFT, there is a need to have some pre-existing knowledge, and the following is the table of contents.
Where i is the imaginary unit i.e., that is, the imaginary unit
Complex form:, where i is the imaginary unit
Complex multiplication: for the two complex sums, then
Because of Euler’s formula (see Equation 1) makes then the complex number
Where is the radius of the circle of the complex plane in which this complex number is located, and is the amplitude angle of this complex number in the complex plane. Then the two complex numbers are, that is,
According to equation 4, the multiplication of two complex numbers can be regarded as the addition of the magnitude angle and the multiplication of the modulus.
Unit Roots: for complex numbers satisfying the equation, we call them nth unit roots. Since according to the multiplication of complex numbers, we know that: multiplication of complex numbers is the addition of magnitudes and the multiplication of moduli. Then for each unit root, the modulus length is 1 and n times the magnitude angle is 0. i.e. (easy to obtain).
To ensure that n times the magnitude angle is always 0, due to this property, we can express the unit root as, where.
We find that n times is always 0 regardless of the value of k.
Recall that, then the n times unit root can be expressed as
Expression of coefficients of a polynomial: given a polynomial
Where is the vector of coefficients
Expression of point values of a polynomial: given a polynomial such as Eqn. (5), we discretize it and set it to take (here why is n) +1 term, will be talked about in section 4) different values from each other, which can be obtained by substituting them, then it is a point-valued expression on.
Multiplication of polynomial coefficient expressions: Given two polynomials
The multiplication of polynomial coefficient expressions is
where:
With Eq. (9), the complexity of the computation is
Multiplication of polynomial point-value expressions: Given two polynomials such as Eqs. (6) and (7), the point-value expressions of their on point-valued expressions are:
Then the multiplication of polynomial point-valued expressions is
See, when we know the complexity in which we can find the point-valued expression of the resultant polynomial.
For the multiplication of a polynomial, according to the above prior knowledge, we can learn: the way to reduce the complexity of polynomial multiplication is to transform the common polynomial coefficients expression into a polynomial point-valued expression, after doing the multiplication of point-valued expression, and then finally transform the point-valued expression into coefficients expression, you can complete the polynomial multiplication.
So the problem becomes:
1. How to transform polynomial coefficient expression into polynomial point-value expression
2. How to transform polynomial point-value expression into polynomial coefficient expression
This leads to the discrete Fourier Transformation (DFT) and the inverse discrete DiscreteFourierTransformation (DFT) and InverseDiscreteFourierTransformation (IDFT)
One way to discretize a polynomial is to substitute values into the polynomial and find the point values in turn. Obviously, the complexity of this approach is yes, which conflicts with our idea of reducing complexity.
So we introduced the Cooley-Tukey algorithm, a classic algorithm for FFT, to reduce the complexity of discretization, which is based on a partitioning strategy.
Given a polynomial
we split it into two polynomials with the same number of terms according to the parity terms (add the polynomials to, with the complementary terms coefficients of 0. PS: why terms, will be mentioned later):
Evidently:
Before going further: we need to return the unit root of the point of knowledge. According to the expression n times unit root, we can obtain an equation
We substitute it into the equation (14), (15) can be obtained:
that is,
So far we find that the original need to substitution of the value of the equation to times reduced to times, and in order to recycle it down, then we only need to recursion times to be able to in the complexity of the That is, we get a point-value pair with an acceptable complexity.
For a more rigorous proof, the following procedure is for those who still have questions
Since Eq. (16) can be obtained
then
The reason for the direct replacement of the sum in the summation is that after squaring, the minus sign is canceled out.
The complexity formula is then
The above is the idea of Cooley-Tukey discrete Fourier transform DFT.
After the DFT, we convert the coefficient expression of the polynomial to the point-valued expression of the polynomial. After completing the multiplication operation, we need to convert the point-valued expression of the polynomial to the coefficient expression of the polynomial in order to obtain the coefficient transformation. The method we use in this case is the inverse discrete Fourier transform IDFT, which is the inverse of the DFT.
The process of solving the IDFT is actually a problem of solving a linear equation, given a linear equation as follows:
The matrix form is as follows:
Assuming that the above matrix is the same as the one above, for the values of the matrices
Set the result of multiplying the two matrices to
Then,
Then,
Then,
At that time,
(where the above equation holds because it is a subunit root, and because the subunit root has a power of 1)
So
then
then the IDFT is the final coefficients that will be obtained by doing the DFT of the result one more time.
In the specific implementation of the FFT, also need to consider for each recursion we how to make a reasonable division. So here the bitreverse algorithm is introduced, also called the butterfly transform.
Through this method of division, we can also summarize another law, for a number of numbers in any position of the number, assuming that the original position for the binary inversion of the function for the last position of its location (the first position of 0), which is a bit binary.
Example:
For a number of figures, then the reversal of its bit binary, then its final position for the first position (ps: the figure does not continue to divide all the figures into one per group, the reader can divide the test)
Here you can add a write-up: if we will be the original array is defined as the reversed array is defined as, then.
And because if it is even, then, then for, but considering that if it is odd, then, then for which is the smallest value satisfied.
The synthesized write-up can be written as
With this write-up, we can directly write the result after DFT division of all numbers.
Reference: from polynomial multiplication to fast Fourier