Power spectrum estimation is widely used in the field of analyzing the frequency components of stationary ergodic random signals, and has been successfully applied to practical projects such as radar signal processing and fault diagnosis. This article gives several types of classical power spectrum estimation methods, and analyzes the performance of the classical power spectrum estimation method through experimental simulation in Matlab. Finally, it explains the limitations of the classical power spectrum estimation method and the reasons for this limitation.
1 Introduction
Given a standard sinusoidal signal, we can analyze its frequency components through the Fourier transform. However, in actual engineering applications, due to the existence of various interferences and noises, the signals we obtain are often not ideal. As shown in Figure 1-1, the signal is uncertain, the amplitude cannot be predicted, and is non-periodic, but it often obeys a certain The statistical characteristics of this signal are called random signals. It should be noted that the random signals mentioned in this article refer to stationary ergodic random signals, and the analysis method of non-stationary random signals [1] will not be discussed in this article.
Figure 1-1 A time domain form of a random signal
For the random signal in Figure 1-1, we can analyze its frequency components through the power spectrum. Figure 1-2 shows the power spectrum of the random signal in Figure 1-1. In the actual process, we can only obtain some discrete data points of the random signal (assumed to be N). This article will discuss how to use these N data points to obtain an “inexact” power spectrum to estimate the power of the real random signal. The power spectrum is estimated and how to achieve better estimation is discussed, that is, several classic power spectrum estimation methods will be described in the next chapter.
Figure 1-2 The random signal power spectrum shown in the above figure
2. Classical power spectrum estimation method
In the previous chapter, we already knew that the power spectrum estimation method uses the N data points that have been obtained to obtain an “inexact” power spectrum to estimate the power spectrum of the real random signal, so before giving the specific method, How to evaluate the quality of this “inaccurate” power spectrum we obtained?
There are many standards for evaluating the performance of power spectrum. This article only gives the two most influential standards: resolution and variance. The resolution is the smallest adjacent frequency component that can be distinguished on the power spectrum. The higher the resolution, the clearer the frequency components of the signal we observe; the size of the variance reflects the volatility of the power spectrum. If the variance is too large, the volatility of the power spectrum will If it is large, it is easy to cause useful frequency components to be drowned by noise. Therefore, for the “inaccurate” power spectrum we hope to obtain, the higher the resolution, the better, and the smaller the variance, the better.
At the same time, we give the concepts of consistent estimation and non-uniform estimation learned in probability theory and mathematical statistics, assuming that the power spectrum of the real signal is
, the estimated “inaccurate” power spectrum is
, if it satisfies formula (2-1), it is called a consistent estimate.
(2-1)
2.1 Periodogram method
2.1.1 Principle of periodogram method
We know N discrete data points
, perform Fourier transform on these data points, and obtain equation (2-2):
(2-2)
Then take the square of the modulus of equation (2-2) and divide it by N to obtain an “inexact” spectrum, as shown in equation (2-3). This is the principle of the periodogram method.
(2-3)
2.1.2 Periodogram method performance (Matlab simulation)
In the previous section we have given the principle of the periodogram method. This section will use Matlab simulation to show the impact of the number of data points N on the performance of the power spectrum. As mentioned above, the resolution and variance of the resulting power spectrum will be analyzed.
We constructed a random signal in Matlab by superposing three sine functions and white noise. Its mathematical form is as formula (2-4), where the frequency
,
,
, amplitude
,
,
, phase
to be independent of each other
The above obeys uniformly distributed random phases,
It is real-valued Gaussian white noise with a mean of 0 and a variance of 1, and the sampling frequency is 1000. The time domain form of the signal is shown in Figure 2-1.
(2-4)
Figure 2-1 Random signals used in the experiment
When the number of data points N is 128, 256, 512 and 1024 respectively, the obtained power spectra are shown in Figure 2-2, Figure 2-3, Figure 2-4 and Figure 2-5 respectively. The resolution can be seen intuitively through the power spectrum graph, and the value of the variance is given in Table 2-1.
Figure 2-2 Power spectrum obtained by periodogram method when N=128
Figure 2-3 The power spectrum obtained by the periodogram method when N=256
Figure 2-4 Power spectrum obtained by periodogram method when N=512
Figure 2-5 Power spectrum obtained by periodogram method when N=1024
Table 2-1 Variance value of power spectrum obtained with different N values
|
N |
128 |
256 |
512 |
1024 |
|
variance |
92.7108 |
130.9109 |
160.9187 |
483.5894 |
By comparing the above experimental results, we can easily find that the resolution and variance of the power spectrum obtained by the periodogram method increase as the number of data points N increases.
2.1.3 Average period diagram method
The power spectrum obtained by the periodogram method is inconsistent with the “high resolution and small variance” we expect. In order to further reduce the variance, the N observation sample data points are
It is divided into L segments, and the data length of each segment is M. The periodogram power spectrum is estimated for each segment of data, and then the average value is obtained. This method is called the average periodogram method.
So how would this approach improve variance? We give proof:
(2-5)
in:
(2-6)
From equation (2-5) we can see that the average period diagram method changes the original variance into the original
, L is the number of segments.
2.1.4 Average periodogram method performance (Matlab simulation)
When the number of data points N is 1024 and the number of segments are 8, 4, and 2 respectively, the power spectra obtained by the average periodogram method are shown in Figure 2-6, Figure 2-7, and Figure 2-8 respectively. The resolution can be seen intuitively through the power spectrum graph, and the value of the variance is given in Table 2-2.
Figure 2-6 Power spectrum obtained by the average periodogram method when L=8
Figure 2-7 Power spectrum obtained by the average periodogram method when L=4
Figure 2-8 Power spectrum obtained by the average periodogram method when L=2
Table 2-2 Variance value of power spectrum obtained with different L values
|
L |
8 |
4 |
2 |
1 |
|
variance |
96.3756 |
190.9647 |
400.6464 |
483.5894 |
When L=1, the average periodogram method degenerates into the periodogram method. By comparing the above experimental results, we can easily find that the power spectrum obtained by the average periodogram method becomes larger as the number of segments L becomes larger, the variance becomes smaller, but the resolution becomes smaller.
When the number of observation sample sequence data N is fixed, to reduce the variance, the number of segments L needs to be increased. When N is not large, the segment length M takes a small value, and the power spectrum resolution is reduced to a lower level. If the number of segments L is fixed, increasing the resolution requires a segment length M, and a longer detection data sequence needs to be collected. In practice, it is precisely the length of the detection sample sequence that is insufficient.
2.1.5 Modified average periodogram method
As mentioned in the previous section, the length of the detection sample sequence is limited in practice. For the existing data length N, if more segments can be divided, a smaller variance will be obtained. Allow overlap between data segments to get more segments. The simplest selection of the overlap length between segments is half of the segment length M. It can be seen from equation (2-5) that more segments can further reduce the variance.
The process of data truncation is equivalent to adding a rectangular window to the data. The larger side lobes of the rectangular window will cause “spectrum leakage”. The window functions we adopt when segmenting are more diverse (triangular windows, Hamming windows, etc.) to reduce the impact of truncated data (adding rectangular windows) window functions [2]
2.1.6 Modified average periodogram method performance (Matlab simulation)
Using the modified average periodogram method, the power spectrum obtained by using the rectangular window, Blackman window and Hamming window respectively is shown in Figure 2-9.
Figure 2-9 Power spectrum obtained by modified average periodogram method with different window functions
It can be found that the rectangular window has the highest resolution, but also the largest variance. This is because the rectangular window spectrum has the narrowest main lobe, the highest resolution, and high side lobes, resulting in the most serious spectrum leakage and the largest variance. MATLAB code classic power spectrum estimation code
2.1.7 Summary
The power spectrum obtained by the periodogram method has a larger resolution and a larger variance as the number of sample points increases; the average periodogram method further improves the variance at the expense of resolution; the modified average periodogram method allows the overlap of segments to further increase the resolution. The number of segments or segments is the same, and the number of sample points in each segment becomes larger. No matter which method is used, there is no complete solution to the contradiction between variance and resolution.
2.2 Correlation power spectrum estimation method-BT method
As we introduced before, to improve the resolution of power spectrum estimation, the length N of the data sequence must be increased, but with longer data sequences, the randomness caused by noise is more fully reflected – larger variance. In fact, when N is infinite, the variance is a non-zero constant. That is, the periodogram method cannot achieve consistent estimation of the power spectrum. The correlation power spectrum estimation method described in this section (hereinafter referred to as the BT method) is a consistent estimation.
2.2.1 Principle of BT method
Wiener Zinchin’s theorem states that the correlation function of a random signal and its power spectrum are a pair of Fourier transforms. The BT method is based on this principle. First, the observation data
The autocorrelation function is estimated, and then the Fourier transform of the autocorrelation function is obtained, and this transformation is used as an estimate of the power spectrum, which is also called the indirect method. The BT method requires that signals other than the signal length N be zero, which also causes the limitations of the BT method.
autocorrelation function
The definition is as shown in Equation (2-7), and the obtained power spectrum is expressed as
, then the BT method can be expressed as formula (2-8).
(2-7)
(2-8)
2.2.2 Performance of BT method (Matlab simulation)
The power spectra obtained by the BT method with the number of data points N being 128, 256, 512 and 1024 are shown in Figure 2-10, Figure 2-11, Figure 2-12 and Figure 2-13.
Figure 2-10 When N=128, the power spectrum obtained by the BT method
Figure 2-11 When N=256, the power spectrum obtained by the BT method
Figure 2-12 When N=512, the power spectrum obtained by the BT method
Figure 2-13 When N=1024, the power spectrum obtained by the BT method
From the above experiments, it can be found that when M increases with the increase of N, the resolution increases and the variance becomes larger. The BT method still does not solve the contradiction between resolution and variance, but when N is infinite, the variance of the power spectrum obtained by the BT method will tend to zero, which is a consistent estimate [2].
2.2.3 The relationship between periodogram method and BT method
related functions
It can be written in the convolution form of equation (2-9):
(2-9)
Let the sequence
The Fourier transform of
, then when M=N-1, the estimation of the power spectrum can be expressed in the form of Equation (2-10). It can be seen that the periodogram method can be regarded as a special case of the BT method when M=N-1.
(2-10)
Conclusion
Through Matlab simulation, this article takes a specific random signal as an example to briefly introduce the basic principles of the periodogram method, average periodogram method, modified average periodogram method and BT method, and analyzes the performance of these methods. It can be seen that neither the periodogram method and its improved algorithm nor the BT method fundamentally solve the contradiction between resolution and variance. Classical power spectrum estimation uses Fourier transform to estimate the power spectrum. However, our previous analysis of random signals does not meet the conditions of Fourier transform, so the classical power spectrum estimation method has to intercept finite-length data points from infinite-length data points and add restrictions. Conditions (the periodogram method actually assumes that the data outside N points are repeated periodically, and the BT method assumes that the data outside N points is zero) to “force” the Fourier transform, which is also the reason for its limitations.
Reference materials
[1] Zhu Zhe, Zhong Hongwei. Research on non-stationary random signal analysis and processing methods [J] Journal of Anhui Institute of Electronic Information Technology 2008.6:28-28
[2] Huangfu Kan. Modern digital signal processing[M]. Electronic Industry Press
Power spectrum estimation
Part of the matlab program code:
Periodogram method: (by classmate Song)
1 Fs=1000;
2 f1=50;
3 f2=125;
4 f3=135;
5 N=128;
6 Nfft=N;
7 n=0:N-1;
8 t=n/Fs;% time series adopted
9 xn=cos(2*pi*f1*t) + 1.5*cos(2*pi*f2*t) + cos(2*pi*f3*t) + 1.5*randn(size(n));
10 figures;
11 plot(n,xn);grid on;title('Time domain signal');
12 P=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier amplitude spectrum squared average value and converted to dB
13 f=(0:length(P)-1)*Fs/length(P);% gives the frequency sequence
14 figures;
15 plot(f(1:N/2),P(1:N/2));grid on;title('Power spectrum (dB chart)');ylabel('Power spectrum/dB');
16 xlabel('frequency/Hz');
17 Pxx=abs(fft(xn,Nfft).^2)/N;%Fourier amplitude spectrum squared average value and converted to dB
18 f=(0:length(Pxx)-1)*Fs/length(Pxx);% gives the frequency sequence
19 figures;
20 plot(f(1:N/2),Pxx(1:N/2));% draw the power spectrum curve
21 xlabel('frequency/Hz');ylabel('power spectrum');
22 title('Periodogram N=128');grid on;
23 std(Pxx)^2
twenty four
25 N=256;
26 Nfft=N;
27 n=0:N-1;
28 t=n/Fs;% time series adopted
29 xn=cos(2*pi*f1*t) + 1.5*cos(2*pi*f2*t) + cos(2*pi*f3*t) + 1.5*randn(size(n));
30 figures;
31 plot(n,xn);grid on;title('Time domain signal');
32 P=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier amplitude spectrum squared average value and converted to dB
33 f=(0:length(P)-1)*Fs/length(P);% gives the frequency sequence
34 figures;
35 plot(f(1:N/2),P(1:N/2));grid on;title('Power spectrum (dB chart)');ylabel('Power spectrum/dB');
36 xlabel('frequency/Hz');
37 Pxx=abs(fft(xn,Nfft).^2)/N;%Fourier amplitude spectrum squared average value and converted to dB
38 f=(0:length(Pxx)-1)*Fs/length(Pxx);% gives the frequency sequence
39 figures;
40 plot(f(1:N/2),Pxx(1:N/2));% draw the power spectrum curve
41 xlabel('frequency/Hz');ylabel('power spectrum');
42 title('Periodogram N=256');grid on;
43 std(Pxx)^2
44
45 N=512;
46 Nfft=N;
47 n=0:N-1;
48 t=n/Fs;% time series adopted
49 xn=cos(2*pi*f1*t) + 1.5*cos(2*pi*f2*t) + cos(2*pi*f3*t) + 1.5*randn(size(n));
50 figures;
51 plot(n,xn);grid on;title('Time domain signal');
52 P=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier amplitude spectrum squared average value and converted to dB
53 f=(0:length(P)-1)*Fs/length(P);% gives the frequency sequence
54 figures;
55 plot(f(1:N/2),P(1:N/2));grid on;title('Power spectrum (dB chart)');ylabel('Power spectrum/dB');
56 xlabel('frequency/Hz');
57 Pxx=abs(fft(xn,Nfft).^2)/N;%Fourier amplitude spectrum squared average value and converted to dB
58 f=(0:length(Pxx)-1)*Fs/length(Pxx);% gives the frequency sequence
59 figures;
60 plot(f(1:N/2),Pxx(1:N/2));% draw the power spectrum curve
61 xlabel('frequency/Hz');ylabel('power spectrum');
62 title('Periodogram N=512');grid on;
63 std(Pxx)^2
64
65 N=1024;
66 Nfft=N;
67 n=0:N-1;
68 t=n/Fs;% time series adopted
69 xn=cos(2*pi*f1*t) + 1.5*cos(2*pi*f2*t) + cos(2*pi*f3*t) + 1.5*randn(size(n));
70 figures;
71 plot(n(1:1000),xn(1:1000));grid on;title('Time domain signal');
72 P=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier amplitude spectrum squared average value and converted to dB
73 f=(0:length(P)-1)*Fs/length(P);% gives the frequency sequence
74 figure;
75 plot(f(1:N/2),P(1:N/2));grid on;title('Power spectrum (dB chart)');ylabel('Power spectrum/dB');
76 xlabel('frequency/Hz');
77 Pxx=abs(fft(xn,Nfft).^2)/N;%Fourier amplitude spectrum squared average value and converted to dB
78 f=(0:length(Pxx)-1)*Fs/length(Pxx);% gives the frequency sequence
79 figures;
80 plot(f(1:N/2),Pxx(1:N/2));% draw the power spectrum curve
81 xlabel('frequency/Hz');ylabel('power spectrum');
82 title('Periodogram N=1024');grid on;
83 std(Pxx)^2
Average period chart method
1 clear;
2Fs=1000;
3 f1=50;
4 f2=125;
5 f3=135;
6 n=0:1/Fs:1;
7 xn=cos(2*pi*f1*n) + 1.5*cos(2*pi*f2*n) + cos(2*pi*f3*n) + 1.5*randn(size(n));
8
9 N=1024;Nsec=512;% data length and data length used by FFT
10 Pxx1=abs(fft(xn(1:512),Nsec).^2)/Nsec; %The first power spectrum
11 Pxx2=abs(fft(xn(513:1000),Nsec).^2)/Nsec;%Second power spectrum
12 Pxx=10*log10((Pxx1 + Pxx2)/2);%Fourier amplitude spectrum squared average value and converted to dB
13 (std((Pxx1 + Pxx2)/2))^2
14 f=(0:length(Pxx)-1)*Fs/length(Pxx);% gives the frequency sequence
15 figures;
16 plot(f(1:Nsec/2),Pxx(1:Nsec/2));% draw the power spectrum curve
17 xlabel('frequency/Hz');ylabel('power spectrum/dB');
18 title('N=2*512');grid on;
19
20 N=1024;Nsec=256;% data length and data length used by FFT
21 Pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec; %The first power spectrum
22 Pxx2=abs(fft(xn(257:512),Nsec).^2)/Nsec;%Second power spectrum
23 Pxx3=abs(fft(xn(513:768),Nsec).^2)/Nsec;%The third power spectrum
24 Pxx4=abs(fft(xn(769:1000),Nsec).^2)/Nsec;%The fourth power spectrum
25 Pxx=10*log10((Pxx1 + Pxx2 + Pxx3 + Pxx4)/4);%Fourier amplitude spectrum squared average value and converted to dB
26 std((Pxx1 + Pxx2 + Pxx3 + Pxx4)/4)^2
27 f=(0:length(Pxx)-1)*Fs/length(Pxx);% gives the frequency sequence
28 figures;
29 plot(f(1:Nsec/2),Pxx(1:Nsec/2));% draw the power spectrum curve
30 xlabel('frequency/Hz');ylabel('power spectrum/dB');
31 title('N=4*256');grid on;
32
33 N=1024;Nsec=128;% data length and data length used by FFT
34 Pxx1=abs(fft(xn(1:128),Nsec).^2)/Nsec; %The first power spectrum
35 Pxx2=abs(fft(xn(129:256),Nsec).^2)/Nsec;%Second power spectrum
36 Pxx3=abs(fft(xn(257:384),Nsec).^2)/Nsec;%The third power spectrum
37 Pxx4=abs(fft(xn(385:512),Nsec).^2)/Nsec;%The fourth power spectrum
38 Pxx5=abs(fft(xn(513:640),Nsec).^2)/Nsec; %The fifth section power spectrum
39 Pxx6=abs(fft(xn(641:768),Nsec).^2)/Nsec;%Sixth section power spectrum
40 Pxx7=abs(fft(xn(769:896),Nsec).^2)/Nsec;%Seventh section power spectrum
41 Pxx8=abs(fft(xn(897:1000),Nsec).^2)/Nsec;%Eighth section power spectrum
42 Pxx=10*log10((Pxx1 + Pxx2 + Pxx3 + Pxx4 + Pxx5 + Pxx6 + Pxx7 + Pxx8)/8);%Fourier amplitude spectrum squared average value and converted to dB
43 std((Pxx1 + Pxx2 + Pxx3 + Pxx4 + Pxx5 + Pxx6 + Pxx7 + Pxx8)/8)^2
44
45 f=(0:length(Pxx)-1)*Fs/length(Pxx);% gives the frequency sequence
46 figures;
47 plot(f(1:Nsec/2),Pxx(1:Nsec/2));% draw the power spectrum curve
48 xlabel('frequency/Hz');ylabel('power spectrum/dB');
49 title('N=8*128');grid on;
Modified average period chart method (by classmate Zhou)
1 clear;% Find the power spectrum and variance of 1024 points
2Fs=1000;
3 f1=50;
4 f2=125;
5 f3=135;
6 n=0:1/Fs:1;
7 xn=cos(2*pi*f1*n) + 1.5*cos(2*pi*f2*n) + cos(2*pi*f3*n) + 1.5*randn(size(n));
8 M=512;N=1024;%The total number of data points is 1024, and the length of each segment is M
9 window=hamming(M);
10 Pxxtotal=0;
11 L=N/(M/2)-1;
12 for i=1:1:(L-1)
13 m=abs(fft(window'.*(xn((M/2 + M/2*i-M + 1):(M/2 + M/2*i))),M).^2)/M ;
14 Pxxtotal=Pxxtotal + m;
15 end
16 window=hamming(Fs-(N-M + 1) + 1);
17 mend=abs(fft(window'.*(xn((N-M + 1):Fs)),M).^2)/M;
18 Pxxtotal=(Pxxtotal + mend)/L;
19 Pxx=10*log10((Pxxtotal));%Fourier amplitude spectrum squared average value and converted to dB
20 f=(0:length(Pxx)-1)*Fs/length(Pxx);% gives the frequency sequence
21 w1=fft(window,N);
22 w10=abs(fftshift(w1));
twenty three
24 plot(f(1:M/2),Pxx(1:M/2));% draw the power spectrum curve
25 xlabel('frequency/Hz');ylabel('power spectrum/dB');
26 grid on;
27 B=var(Pxxtotal)
BT method (by Yang Yang)
1 Fs=1000;
2 n=0:1/Fs:1;
3 x=cos(2*pi*50*n) + 1.5*cos(2*pi*125*n) + cos(2*pi*135*n) + 1.5*randn(size(n));
4nfft=1024;
5 ncxk=3*nfft/4;
6 cxn=xcorr(x,'unbiased');
7 CXk=fft(cxn,ncxk);
8 Pxx=abs(CXk);
9 index=0:round(ncxk/2-1);
10 k=index*Fs/ncxk;
11 C=Pxx(index + 1);
12 P=(Pxx(index + 1));
13 plot(k,P);grid on
14var(C)
15 title('BT method power spectrum estimation N=1024');
16 xlabel('frequency Hz');
17 ylabel('power');