那位朋友能帮我解释看些这段程序啊

[shuijnglei0125]

那位MATLAB高手能帮我看下下面的程序啊、

帮忙说下整体程序实现的功能

最好能给没句话解释下



%求信号x的功率谱
%时域信号--(FFT变换)-->幅度谱--(平方)-->功率谱 直接法
%时域信号--->相关函数--(FFT变换)-->功率谱--(除以频率分辨率)-->功率谱密度 间接法
clear all
%a函数及其频谱
m=-20:0.2:20;
index1=find(m<0);
a = 10*sin(m(index1))./(m(index1));
a=[a 10];
index2 = find(m > 0);
a = [a 10*sin(m(index2))./(m(index2))];
figure(1);subplot(411);stem(m,a);title('原函数a');
M=length(m);
p=0:M-1;q=0:M-1;
A=a*exp(-1i*2*pi/M).^(p'*q);
subplot(412);plot(fftshift(abs(A)));title('原函数a的幅度谱');
%调制函数b
t=0:1/100:2;
b=100*sin(2*pi*20*t);
f=linspace(0,100,128);
subplot(413);plot(t,b);title('其中的一个调制函数b=50sin(2pi20t)');
B=fft(b,128);
subplot(414);plot(f,abs(B));title('其中的一个调制函数b的幅度谱');
%调制函数c
c=80*sin(2*pi*40*t);
ba=a.*b;ca=a.*c;
x=ba+ca+randn(size(ba));
X=abs(fft(x,128));
nfft=512;
index=0:nfft-1;
k=index*100/nfft;%100采样频率
[r tao]=xcorr(x,'unbiased');
subplot(4,1,2);plot(tao,r);title('x的自相关函数');
Pxx=abs(fft(r,nfft));
Sf=10*log10(Pxx(index+1)); %功率谱密度函数
n0=length(X);
n1=length(Sf);
% 多尺度
m=-4:-1;
for i=1:4
delta(i)=2.^m(i);
end
N=16; % 滤波器长度（需要调整，必须是偶数）
A=-1/sqrt(2*pi); % 幅度
% 构造高斯函数
for index_x=1:N;
x=index_x-(N+1)/2;
W1(index_x)=A*(x/delta(1)^2).*exp(-(x.*x)/(2*delta(1)^2));
W2(index_x)=A*(x/delta(2)^2).*exp(-(x.*x)/(2*delta(2)^2));
W3(index_x)=A*(x/delta(3)^2).*exp(-(x.*x)/(2*delta(3)^2));
W4(index_x)=A*(x/delta(4)^2).*exp(-(x.*x)/(2*delta(4)^2));
end;
W1=W1./norm(W1);W2=W2./norm(W2);W3=W3./norm(W3);W4=W4./norm(W4); % 能量归一化
%对信号的幅度谱做奇异性检测
g0=conv(X,W1).*conv(X,W2).*conv(X,W3).*conv(X,W4);
g0=wkeep(g0,n0); % 保持信号长度
m0=abs(g0); % 取模
[V0,P0]=findpeaks(m0); % 模极大值
F0=f(P0);
figure(2);
subplot(2,1,1);plot(f,X);
xlabel('频率/Hz');ylabel('信号幅值');title('原始信号');
subplot(2,1,2);
hold on;
plot(f,m0);
plot(F0,V0,'r*');
title(['小波模最大值对应的频率是',num2str(F0)]);
xlabel('频率/Hz');ylabel('小波变换后的模值');
%对信号的功率谱密度做奇异性检测
g1=conv(Sf,W1).*conv(Sf,W2).*conv(Sf,W3).*conv(Sf,W4);
g1=wkeep(g1,n1);
m1=abs(g1);
[V1,P1]=findpeaks(m1);
F1=k(P1);
figure(3);
subplot(2,1,1);plot(k,Sf);
xlabel('频率/Hz');ylabel('信号幅值');title('原始信号');
subplot(2,1,2);
hold on;
plot(k,m1);
plot(F1,V1,'r*');
title(['小波模最大值对应的频率是',num2str(F1)]);
xlabel('频率/Hz');ylabel('小波变换后的模值');