gonglvpu.zip_periodogram window_welch 窗函数_信号平均_平均功率法_窗函数

%% 产生信号 fs = 800; % 采样频率 t = 0:1/fs:0.5; x = cos(2*pi*t*200) + 1.5*cos(2*pi*t*250) + normrnd(0,1,(size(t))); % 随机信号 两个正弦信号和一个高斯白噪声 figure(1); subplot(2,1,1); plot(t,x), xlabel('t'), ylabel('x'); title('输入信号（时域）');grid on; %画出时域图形 %% 1、 BT法 win1 = rectwin(length(x)); % 得到点数N等于x1长度的矩形窗 xn = x'.*win1; % x1与窗函数卷积 Rx = xcorr(xn); % 求xn的自相关函数Rx N = length(Rx); %自相关函数的长度 Sw = abs(fft(Rx, N)); %对Rx做N点的FFT得到功率谱密度Sw n = 0:N-1; k = 0:N-1; %画出Rx、Sw及归一化Sw的图像 figure(1); subplot(2,1,2); plot(n*fs/N, Rx(n+1),'r'); title('BT法 自相关函数Rx'); xlabel('t2-t1');ylabel('Rx'),grid on; hold on; figure(2); subplot(2,1,1); plot(k*fs/N, Sw(k+1),'g'); title('BT法 功率谱估计'); xlabel('f(Hz)');ylabel('功率谱密度'),grid on; hold on; figure(2); subplot(2,1,2); plot(k*fs/N, 10*log10(Sw(k+1)/max(Sw)),'y'); title('BT法 归一化功率谱估计'); xlabel('f(Hz)');ylabel('归一化功率谱密度(dB)'),grid on; hold on; %% 2、周期图法 win_2 = rectwin(length(x)); % 得到点数N等于x1长度401的矩形窗 xn = x'.*win_2; index = nextpow2(length(x)); %最靠近且比length(x)大的2次幂的幂数 N = pow2(index); %求最靠近且比length(x)大的2的幂 Xw = fft(xn, N); % 对Xn做N点FFT Sw = Xw.*conj(Xw)/N; % 计算功率谱密度 conj（）：求共轭 k = 0:N-1; n = 0:N-1; figure(3); subplot(2,1,1); plot(n*fs/N, Sw(n+1)); title('周期图法 功率谱估计'); xlabel('f(Hz)');ylabel('功率谱密度'),grid on; hold on; subplot(2,1,2); plot(k*fs/N, 10*log10(Sw(k+1)/max(Sw)),'g'); title('周期图法 归一化功率谱估计'); xlabel('f(Hz)');ylabel('归一化功率谱密度(dB)'),grid on; hold on; %% 3、窗函数法 win3 = rectwin(length(x)); % 得到点数N等于x1长度401的矩形窗 xn = x'.*win3; % x1与窗函数卷积 Rx = xcorr(xn); % 求xn的自相关函数Rx 801个点 N = length(Rx); %自相关函数的长度 win_4 = bartlett(length(N)); Rxx = Rx*win_4; Sw = abs(fft(Rxx, N)); %对Rx做N点的FFT得到功率谱密度Sw n = 0:N-1; k = 0:N-1; %画出Rx、Sw及归一化Sw的图像 figure(4); subplot(2,1,1); plot(n*fs/N, Sw(n+1)); title('窗函数法 功率谱估计'); xlabel('f(Hz)');ylabel('功率谱密度'),grid on; hold on; subplot(2,1,2); plot(k*fs/N, 10*log10(Sw(k+1)/max(Sw)),'y'); title('窗函数法 归一化功率谱估计'); xlabel('f(Hz)');ylabel('归一化功率谱密度(dB)'),grid on; hold on; %% 4、平均法 segment4 = 50; % 每一段的长度 num = 8; % 总共分成的段数 Sww = 0; xm = zeros(1,50); for i = 1 : num window_c = rectwin(50); for j = 1 : 50 xm(j) = x(j+(i - 1)*50); end N = 64; %最靠近且比segment大的2的幂数 即64 Xw = fft(xm, N); % 对Xn做N点FFT Sw = Xw.*conj(Xw)/N; % 计算功率谱密度 Sww = Sww + Sw; end Sww = Sww/10; k = 0:N-1; n = 0:N-1; figure(5); subplot(2,1,1); plot(n*fs/N, Sww(n+1)); title('平均法 功率谱估计'); xlabel('f(Hz)');ylabel('功率谱密度'),grid on; hold on; subplot(2,1,2); plot(k*fs/N, 10*log10(Sww(k+1)/max(Sww)),'g'); title('平均法 归一化功率谱估计'); xlabel('f(Hz)');ylabel('归一化功率谱密度(dB)'),grid on; hold on; %% 5、Welch法 segment5 = 64; % 每一段的长度 noverlap = 16; % 两端重叠的部分 nfft = pow2(nextpow2(length(x))); % 对每一段Xn做FFT [Sxx,f] = pwelch(x,segment5,noverlap,nfft,fs); % 调用pwelch函数 figure(6); subplot(2,1,1); plot(f,Sxx,'b'); title('Welch法 功率谱估计'); xlabel('f(Hz)');ylabel('功率谱密度'),grid on; hold on; subplot(2,1,2); plot(f,10*log10(Sxx/max(Sxx)),'r'); title('Welch法 功率谱估计'); xlabel('f(Hz)');ylabel('归一化功率谱密度(dB)'),grid on; hold on;