| 几种常用的经典频谱估计方法原理介绍及Matlab算法实现 | 您所在的位置：网站首页 › 频谱分析技术有哪些类型的 › 几种常用的经典频谱估计方法原理介绍及Matlab算法实现 |
几种常用的经典频谱估计方法原理介绍及Matlab算法实现
|
autocorr_signal_noise = xcorr(signal_noise); %求线性自相关
fft_autocorr_signal_noise = fft(autocorr_signal_noise); %求自相关的fft
abs_fft_autocorr_signal_noise = abs(fft_autocorr_signal_noise)/(2*points-1); %求幅值
max_abs_fft_autocorr_signal_noise = max(abs_fft_autocorr_signal_noise);
abs_fft_autocorr_signal_noise = abs_fft_autocorr_signal_noise/max_abs_fft_autocorr_signal_noise;
log10_abs_fft_autocorr_signal_noise = 10*log10(abs_fft_autocorr_signal_noise);
Delta_f1 = (0:2*points-2)*fs/(2*points-1);
figure;
plot(Delta_f1(1:points)/1000,log10_abs_fft_autocorr_signal_noise(1:points));
xlabel('KHz');
ylabel('归一化功率谱P(k)/dB');
string = ['间接法（自相关法BT）信号点数=',num2str(signal_points),',采样点数=',num2str(points),...
',信号频率=',num2str(f0/1000),'KHz,采样率=',num2str(fs/1000),'KHz,信噪比SNR=',...
num2str(SNR),'dB'];
title(string);
Bartlett法的Matlab程序
首先将截断的信号分成等长的不重叠的若干段，然后计算每一段的离散傅里叶变换，并取模值的平方，然后对应点求和。所有段计算并求和完成后，归一化处理，并取以10为底的对数然后乘以10，绘制归一化后的功率谱图。 M = 128; %每段的点数 L = points/M; %段数 Bartlett_Spectrum = zeros(1,points); for i=1:L temp_signal = zeros(1,points); %清零 start_point = (i-1)*M; temp_signal(1:M) = signal_noise((start_point+1):(start_point+M)); %截取第i段信号 fft_temp_signal = fft(temp_signal); squar_fft_temp_signal = fft_temp_signal.*conj(fft_temp_signal); %求当前段的功率谱 Bartlett_Spectrum = Bartlett_Spectrum + squar_fft_temp_signal; %累加 end max_Bartlett_Spectrum = max(Bartlett_Spectrum); Bartlett_Spectrum = Bartlett_Spectrum/max_Bartlett_Spectrum; %归一化 log10_Bartlett_Spectrum = 10*log10(Bartlett_Spectrum); %求对数 Delta_f2 = (0:points-1)*fs/points; figure; plot(Delta_f2(1:points/2)/1000,log10_Bartlett_Spectrum(1:points/2)); xlabel('KHz'); ylabel('归一化功率谱P(k)/dB'); string = ['Bartlett法,信号点数=',num2str(signal_points),',采样点数=',num2str(points),... ',信号频率=',num2str(f0/1000),'KHz,采样率=',num2str(fs/1000),'KHz,信噪比SNR=',... num2str(SNR),'dB,段数L=',num2str(L),',每段点数M=',num2str(M)]; title(string); Welch法的Matlab程序首先将截断的信号分成等长的有重叠的若干段，然后计算每一段的离散傅里叶变换，并取模值的平方，然后对应点求和。所有段计算并求和完成后，归一化处理，并取以10为底的对数然后乘以10，绘制归一化后的功率谱图。 M = 128; %每段的点数 Overlap = M-M/16; %重叠点数 L = (points - Overlap)/(M-Overlap); %段数 Welch_Spectrum = zeros(1,points); for i=1:L temp_signal = zeros(1,points); %清零 start_point = (i-1)*(M-Overlap); temp_signal(1:M) = signal_noise((start_point+1):(start_point+M)); %截取第i段信号 fft_temp_signal = fft(temp_signal); squar_fft_temp_signal = fft_temp_signal.*conj(fft_temp_signal); %求当前段的功率谱 Welch_Spectrum = Welch_Spectrum + squar_fft_temp_signal; %累加 end max_Welch_Spectrum = max(Welch_Spectrum); Welch_Spectrum = Welch_Spectrum/max_Welch_Spectrum; %归一化 log10_Welch_Spectrum = 10*log10(Welch_Spectrum); %求对数 Delta_f3 = (0:points-1)*fs/points; figure; plot(Delta_f3(1:points/2)/1000,log10_Welch_Spectrum(1:points/2)); xlabel('KHz'); ylabel('归一化功率谱P(k)/dB'); string = ['Welch法,信号点数=',num2str(signal_points),',采样点数=',num2str(points),... ',信号频率=',num2str(f0/1000),'KHz,采样率=',num2str(fs/1000),'KHz,信噪比SNR=',... num2str(SNR),'dB,段数L=',num2str(L),',每段点数M=',num2str(M),',重叠点数=',... num2str(Overlap)]; title(string); Matlab运行结果及分析设置相应参数运行上述程序得到如下图： 
图1信号的时域图
图2周期图法得到的功率谱图 图3 BT法得到的功率谱(M=128)
图4 Batlett法得到的功率谱（4段） 图5 Welch法得到的功率谱(128点)
图6 BT法得到的功率谱(M=64) 图7 BT法得到的功率谱(M=32)
图8 Batlett法得到的功率谱（8段） 图9 Welch法得到的功率谱(64点)
图1为生成的仿真信号，一共512个点。 图2为用周期图法对数据直接求出功率谱。由于主瓣宽度B = 2/512 = 0.003906 < (150-140)/750 = 0.013所以f0与f1能够分开。 图3是利用BT法求出的功率谱，M=128,依然能将f0与f1分开；对于与图2其主瓣更宽了。当M=64时，如图6 f0与f1不能完全分开，只是在波形的顶部看出是两个频率的分量。当M=32,如图7 f0与f1就已经无法分开了。 图4是利用Batlett法求出的功率谱，共分4段，每段128个点，可以看到f0与f1能够被分开，相对于图2其更加平滑，但是主瓣的宽度比图2的要宽。当分为8段，每段64个点时，如图8其变的更加平滑但f0与f1已经不能完全分开了。 图5是利用Welch方法求出的平均周期图，每段128个点，重合120个点，谱变得更加的平滑，效果与图4基本一致。若每段64个点，重合60个点，如图9 f0与f1已不能分开，效果与图8基本一致。 |
【本文地址】