修正的周期图法matlab代码

修正的周期图法

clc;
clear all;
% cd D:\Users\ASUS\Desktop\SIGNAL\signal_1
% load DSPINPUTSIG_MUSIC.mat;
% S1=musicsignal';
% S2=load ('DSPINPUTSIG_TEST.txt')';
Fs=1000;
f1=50;
f2=125;
f3=165;

N=256;%序列长度
n=0:N-1;
t=n/Fs;%采用的时间序列

xn=cos(2*pi*f1*t)+1.5*cos(2*pi*f2*t)+cos(2*pi*f3*t)+randn(size(n));
%%
S=xn;
N=length(S);%数据总长度
Nfft=1024;

%U为修正因子
U=1;
L=64;%L为每段数据长度
R=32;%R为一跳长度
num_L=fix(((N-L)/R)+1)%分段数

matrix_x=zeros(num_L,L);

fft_x=zeros(num_L,Nfft);
I=zeros(num_L,Nfft); 
Ip=zeros(1,Nfft);

for r=1:num_L
    start_xn=(r-1)*R+1;
    end_xn=(r-1)*R+L;
%     start_xn=round((r-1)*R)+1;%每段的第一个数字
%     end_xn=round((r-1)*R+L);%每段的最后一个数字
    
    if end_xn>N             %超出总长度
        matrix_x(r,:)=[S(start_xn:N),zeros(1,end_xn-N)];%超出位数补0
        fft_x(r,:)=fft(matrix_x(r,:),Nfft);%fft变换
        I(r,:)=fft_x(r,:).*conj(fft_x(r,:))/(L*U);%每一段的功率谱
        Ip=Ip+I(r,:);%功率谱相加
    else
    matrix_x(r,:)=S(start_xn:end_xn);
    fft_x(r,:)=fft(matrix_x(r,:),Nfft);
    I(r,:)=fft_x(r,:).*conj(fft_x(r,:))/(L*U);
    Ip=Ip+I(r,:);
    end
end

Ip=Ip/num_L;%求均值
var(Ip)
Ip=Ip/max(Ip);%归一化


figure(1);
p=0:2/Nfft:2-2/Nfft;
plot(p,Ip,'LineWidth',.2);

figure(2);
w=-1:2/Nfft:1-2/Nfft;
plot(w,fftshift(Ip),'LineWidth',.2);
xlabel('数字角频率/pi*rad');

title('INPUTSIG-MUSIC的功率谱估计结果（平均周期图法）','FontSize',10);
grid on;

num_L = 7分为七段
var= 38.0873 方差大小

周期图法

clc;
clear all;

Fs=1000;
f1=50;
f2=125;
f3=165;

N=256;%序列长度
Nfft=1024;
n=0:N-1;
t=n/Fs;%采用的时间序列

xn=cos(2*pi*f1*t)+1.5*cos(2*pi*f2*t)+cos(2*pi*f3*t)+randn(size(n));
%%
S=xn;
N=length(S);%数据总长度

U=1;%U为修正因子
L=64;%L为每段数据长度
R=64;%R为一跳长度

x=zeros(fix(N/R),L);
fft_x=zeros(fix(N/R),Nfft);
I=zeros(fix(N/R),Nfft); 
Ip=zeros(1,Nfft);
num_L=fix(N/R)
for r=1:fix(N/R)
    start_xn=round((r-1)*R)+1;%每段的第一个数字
    end_xn=round((r-1)*R+L);%每段的最后一个数字
    if end_xn>N             %超出总长度
        x(r,:)=[S(start_xn:N),zeros(1,end_xn-N)];%超出位数补0
        fft_x(r,:)=fft(x(r,:),Nfft);%fft变换
        I(r,:)=fft_x(r,:).*conj(fft_x(r,:))/(L*U);%每一段的功率谱
        Ip=Ip+I(r,:);%功率谱相加
    else
    x(r,:)=S(start_xn:end_xn);
    fft_x(r,:)=fft(x(r,:),Nfft);
    I(r,:)=fft_x(r,:).*conj(fft_x(r,:))/(L*U);
    Ip=Ip+I(r,:);
    end
end

Ip=Ip/fix(N/R);%求均值
var(Ip)
Ip=Ip/max(Ip);%归一化
w=-1:2/Nfft:1-2/Nfft;

figure(1);
p=0:2/Nfft:2-2/Nfft;
plot(p,Ip,'LineWidth',.2);
xlabel('数字角频率/pi*rad');
title('INPUTSIG-MUSIC的功率谱估计结果（平均周期图法）','FontSize',10);

figure(2);
w=-1:2/Nfft:1-2/Nfft;
plot(w,fftshift(Ip),'LineWidth',.2);
xlabel('数字角频率/pi*rad');
title('INPUTSIG-MUSIC的功率谱估计结果（平均周期图法）','FontSize',10);
grid on;

num_L =

 4

var =

51.3131方差大