Tip: after improvement is a follow-up to before improvement

preface

In the periodic graph method and BT method, the original signal of power spectrum estimation is set as Gaussian white noise. Because BT method must be used in generalized stationary random process, the existing function is only the function of periodic graph method. Therefore, in the first half, I use a single sample to represent the two methods. Before improvement, the matlab programs are written from the original formula, This makes it easier to understand the internal processes of the two methods. Finally, the quality of power spectrum estimation: deviation and variance will also be given.

In the improved method, we directly use Pwelch function to realize the piecewise periodic graph method; Since there is no direct formula for the segmented BT method (perhaps because the BT method can only be used in the generalized stationary process, or the effect of the BT method is not as good as the periodic method), on the basis of a single sample, the segmented BT method is realized by first processing each segment, then processing each segment with the BT method, and finally averaging. The third method is to balance the resolution and variance by using MTM method.

1, Power spectrum estimation method (before improvement)

1. Periodogram method

2.BT method

3. Relationship between periodogram method and BT method

4. Quality of power spectrum estimation

These spectral estimation methods before improvement have been described in the previous section, and this section introduces the improved methods.

2, Improvement of power spectrum estimation

Segmentation methods are mainly divided into two categories: Bartlett average method and Welch average method. Bartlett average method is a non overlapping segmentation method, while Welch average method is an overlapping segmentation method.
1. The segmentation method based on the periodic graph method can be described by Pwelch function,
Some programs will use psd function, but I help ed it. It is shown as follows. The description is mainly based on pwelch function!

2. The segmentation method based on BT method has no function that can be called directly. This paper analyzes it from the basic formula by combining BT method and two segmentation methods, and also knows the difference between non overlapping segmentation and overlapping segmentation.

In order to compare the various methods, the following procedures are based on x There are 20000 points,
FFT 20000 points, Fs=1; n=0:N-1，And the formula is:
x=cos(0.3*pi*n)+cos(0.32*pi*n)+randn(size(n));

1.Pwelch method (subsection method under periodic graph method)

This method can directly call the function Pwelch. This article calls:
[pxx,f] = pwelch(x,window,noverlap,Nfft,Fs)
The meaning of input parameters is:

• x: Estimated signal

• Window: window function added to smooth the power spectrum. hamming window is used by default
  Window function query

• noverlap: indicates the number of overlaps between two segments

• Nfft: indicates the number of points of fft, nfft > x, and the default value is 256 points

• FS: represents the sampling frequency, and the actual frequency f(Hz)/Fs = normalized frequency w(rad/s)/2pi

Output parameters:

• pxx: power spectrum
• f: Abscissa, range [0 pi*Fs]

If you want to see other forms, help Just a minute.

If boxcar window is used and NOVERLAP=0, the average periodogram of Bartlett method can be obtained.

%% Segmentation method of periodic graph method
clc;close all;clear;
N=20000; Nfft=20000; Fs=1; n=0:N-1;
x=cos(0.3*pi*n)+cos(0.32*pi*n)+randn(size(n));
M=256;                  %%be careful M If the value of is too small, the separation rate will be too small
%% Bartlett average method
window1=boxcar(M); noverlap1=0;        
figure;
subplot(2,1,1);
[Pxx1,f]=pwelch(x,window1,noverlap1,Nfft);
plot(f/pi,Pxx1);grid on;title('Based on periodic graph method Bartlett average method');
%In order to easily see the line spectrum position( pi The previous factor), will f divide pi. 
%If used boxcar Windows, and NOVERLAP=0，You can get Bartlett Average periodogram of method
%% Welch average
window2=hamming(M); noverlap2=M/2;
[Pxx2,f]=pwelch(x,window2,noverlap2,Nfft);
subplot(2,1,2);
plot(f/pi,Pxx2);grid on;title('Based on periodic graph method Bartlett average method');

It is worth noting that the number of points of the window function cannot be set too small, otherwise the resolution may be too small to distinguish the two frequencies. When M=64, the results are as follows:

2. Segmentation method under BT method

Since there is no piecewise function based on BT method, the following is to segment it artificially, and then average the results. In order to do a group of control experiments, a group of non segmented numbers is added (but this number uses one of the segments after segmentation. In fact, when the number of points is unfair, the resolution will not be as large as the number of points, but the variance will be large when there are more points).

%% BT Subsection method of method
clc;clear all;close all;
N=20000;              %%When the number of segments is small, the resolution may not be enough and will increase N Just click
K=200;                %%Number of segments
M=N/K;                %%Sampling points per segment
n=1:N;
x1=cos(0.3*pi*n)+cos(0.32*pi*n)+randn(size(n));
for k=1:K
    x(k,:)=x1(M*k-M+1:M*k);
end
%% Bartlett average method
for i=1:K
    [rx_BT(i,:),lags]=xcorr(x(i,:),'coeff');
    sf_BT(i,:)=fftshift(abs(fft(rx_BT(i,:))));
end
rx_BT_AV=sum(rx_BT)/K;
sf_BT_AV=sum(sf_BT)/K;
lags1=lags/M;       %%Make the horizontal axis in[0 0.5]
figure;
subplot(3,2,1);
plot(lags,rx_BT(1,:));title('Non segmented autocorrelation function');
subplot(3,2,2);
plot(lags1,sf_BT(1,:));title('Non segmented power spectral density');grid on;
axis([0 +inf -inf +inf]);
subplot(3,2,3);
plot(lags,rx_BT_AV);title('be based on BT Legal Bartlett Piecewise average autocorrelation function');
subplot(3,2,4);
plot(lags1,sf_BT_AV);title('be based on BT Legal Bartlett Piecewise average power spectral density');grid on;
axis([0 +inf -inf +inf]);

%% Welch average
for i=1:1:K-1
    [rx_BT1(2*i-1,:),lags]=xcorr(x(i,:),'coeff');
    [rx_BT1(2*i,:),lags]=xcorr([x(i,M/2+1:M),x(i+1,1:M/2)],'coeff');
end
    [rx_BT1(2*K-1,:),lags]=xcorr(x(K,:),'coeff');
    [rx_BT1(2*K,:),lags]=xcorr([x(K,M/2+1:M),x(1,1:M/2)],'coeff');   
for i=1:2*K
        sf_BT1(i,:)=fftshift(abs(fft(rx_BT1(i,:))));
end
rx_BT1_AV=sum(rx_BT1)/(2*K);
sf_BT1_AV=sum(sf_BT1)/(2*K);
subplot(3,2,5);
plot(lags,rx_BT1_AV);title('be based on BT Legal Welch Piecewise average autocorrelation function');
subplot(3,2,6);
plot(lags1,sf_BT1_AV);title('be based on BT Legal Welch Piecewise average power spectral density');grid on;
axis([0 +inf -inf +inf]);

3.MTM method

MTM method is also called multi window method. Instead of using band-pass filters, it uses a set of optimal filters to calculate the estimated value. These optimal FIR filters are obtained from a set of discrete flat sphere like sequence (DPSS, also known as Slepian).
In the function, MTM provides a time bandwidth parameter (the product of time and bandwidth NM), which can be used to balance variance and resolution. The larger the NM, the more power spectrum estimates and the smaller the variance. Each group of data can find the most suitable parameter to balance variance and resolution according to their own needs.

%% MTM method
clc;clear all;close all;
N=20000; Nfft=20000; Fs=1; n=0:N-1;
x=cos(0.3*pi*n)+cos(0.32*pi*n)+randn(size(n));
[pxx1,f]=pmtm(x,4,Nfft,Fs);
figure;
subplot(2,1,1);plot(f,pxx1);title('Multi window method( MTM)NW=4');grid on;
[pxx2,f]=pmtm(x,2,Nfft,Fs);
subplot(2,1,2);plot(f,pxx2);title('Multi window method( MTM)NW=2');grid on;

summary

It can be seen from the power spectrum of the above three methods:
1. The estimation quality based on BT method is the worst (no wonder there is no function to represent the estimation quality based on BT method), and when N is reduced to 2000, two frequencies cannot be identified respectively (it can be copied and modified for verification). The estimation quality of MTM method is the best.
2. The number of segments of Bartlett average method and Welch average method is different, but the effect is not very different, but the number of window functions in pwelch function will affect the resolution.