mypsd.m

现代信号处理里的一些内容

%power spectral density (PSD) estimate
function [Pxx,w]=MyPSD(x,N,L);
if nargin==1
    N=length(x);
    L=1;
elseif nargin==2
    L=1;
end
M=floor(N/L);%保证M为正数；
m=mod(N,L);%N/L的余数，即被舍弃的几个数据
n=[0:M-1];
w=2*pi.*(n-M/2)./M;
X2=[];

subplot(2,2,1);
plot([0:N-1],x(1:N));
title(['原始信号波形(N=',num2str(N),')']);
xlabel('n');
ylabel('x(n)');
grid on;

subplot(2,2,2);
for i=0:L-1
    k=x(i*M+1:(i+1)*M);
    X2=[X2;k];
end
plot(n,X2)
title(['信号分段波形(L=',num2str(L),')']);
xlabel('n');
ylabel('x_{i}(n)');
grid on;

subplot(2,2,3);
Fx=fft(X2,[],2);
Fx=fftshift(Fx,2);
P=(abs(Fx).^2)./M;
%disp('size(w)');
%size(w),
%disp('size(P)');
%size(P)
plot(w,P);
title(['功率谱密度(L=',num2str(L),'),(N=',num2str(N),')']);
xlabel('\varpi');
ylabel('I_{i}(\varpi)');
grid on;

subplot(2,2,4);
Pxx=sum(P,1)./L;
plot(w,Pxx);
title(['功率谱密度(L=',num2str(L),'),(N=',num2str(N),')']);
xlabel('\varpi');
ylabel('P_{xx}(\varpi)');
grid on;