数字信号处理仿真，经典谱估计matlab代码，含注释，简单易懂

包含经典谱估计的多种方法：直接法（周期图法），间接法（自相关法），直接法的改进方法（Bartlett法，Welch法）

clcclear
n=0:255;
fs=1;
x=cos（0.35*pi*n）+cos（0.4*pi*n）+1.5*cos（0.8*pi*n）;
a=0.5;
noise = a*randn（size（x））;                %加高斯白噪声的信号
y = x+ noise;
N=256;

%%%%%%%%%%%%%%%%%%%经典谱估计%%%%%%%%%%%%%%%%%%%%%%%
%%%周期图：对信号做256点的傅里叶变换，取幅值的平方再除以N
F1=abs（fft（yN））;
P_per=1/N*F1.^2;
f=-0.5:1/N:0.5-1/N;
nn=2*pi*f/pi;                    %0:1/fs:（1-1/fs）;
subplot（221）
P_per=10*log10（P_per）;
y2=[P_per（129:256）P_per（1:128）];
plot（nny2）
xlabel（‘频率‘）ylabel（‘功率谱P（dB）‘）title（‘周期图法‘ ）; 

%%%%自相关法（BT法）：对信号求自相关函数再对他做傅里叶变换（平滑）
%****自相关函数的直接计算
rxk=xcorr（y‘unbiased‘）;  
%N=length（rxk）;
pxw21=abs（fft（rxkN））;
subplot（222）
pxw22=10*log10（pxw21）;
y3=[pxw22（129:256）pxw22（1:128）];
plot（nny3）
xlabel（‘频率‘）ylabel（‘功率谱P（dB）‘）title（‘BT法‘ ）; 
%%%%%%%%%%%%Baetlett法
M=32;                                                      %段长
L1=length（y）/M;                                            %段数
Win1=boxcar（M）;                                            %矩形窗
sum_Per=0;
y_ba=[L1M];Y_ba=[L1M];P_per=[L1M];
for i=1:L1
    y_ba（i1:M）=y（1+（i-1）*M:M+（i-1）*M）.*Win1‘;          %加窗
    Y_ba（i1:M）=abs（fft（y_ba（i1:M）M））;                %一段傅里叶变换的幅值
    P_per（i1:M）=1/M*Y_ba（i1:M）.^2;                    %一段的功率谱
end
AA=sum（P_per1）;
AA=10*log10（AA）;
y3=[AA（length（AA）/2+1:length（AA））AA（1:length（AA）/2）];
ff=-0.5:1/M:0.5-1/M;
nn=2*pi*ff/pi;
subplot（223）
plot（nny3）;
xlabel（‘频率‘）ylabel（‘功率谱P（dB）‘）title（[‘ Baetlett法段长M=‘ int2str（M）] ）; 
%%%%%%%%%%%%%%Welch法
M=32;
L2=（N-M/2）/M*2;
Win2=hamming（M）;                  %
kk=Win2.^2;
U=1/M*sum（kk）;                                              %归一化因子
y_we=zeros（L2M）;Y_we=zeros（L2M）;P_per_we=zeros（L2M）;
for i=1:L2
    y_we（i1:M）=y（1+（i-1）*M/2:M+（i-1）*M/2）.*Win2‘;           %加窗
    Y_we（i1:M）=abs（fft（y_we（i1:M）M））;                     %一段傅里叶变换的幅值
    P_per_we（i1:M）=1/M/U*Y_we（i1:M）.^2;                    %一段的功率谱
end
BB=sum（P_per_we1）;
BB=10*log10（BB）;
y4=[BB（length（BB）/2+1:length（BB））BB（1:length（BB）/2）];
ff=-0.5:1/M:0.5-1/M; 
subplot（224）
plot（nny4）;
xlabel（‘频率‘）ylabel（‘功率谱P（dB）‘）title（[‘ Welch法段长M=‘ int2str（M）] ）; 

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        2295  2018-11-23 08:09  puguji.m