微弱信号的检测提取这个程序的提取算法有错误吗,最后得出的波形不大对啊
对微弱信号进行检测提取与分析,最后得出的波形和给出来的波形有的差别很大啊,小波分析消除微弱信号的算法怎么做啊==fs=44100;%设定采样频率
N=1024;%取的样本点
n=0:N-1;
t=n/fs;%采样
x1=wavread('xn1100.wav');%读入混合信号
%******************************************
m1=mean(x1); %求均值
v1=var(x1); %求方差
w1=v1+m1^2; %求均方值
%******************************************
%进行FFT变换并做频谱图
y=fft(x1,N);%进行fft变换
magy=abs(y);%求幅值
angley=angle(y);%求相位
fy=(0:length(y)-1)*fs/length(y);%进行对应的频率转换
figure(1);
plot(fy,magy,'k',fy,angley,'r');%做频谱图
xlabel('频率(Hz)');
ylabel('幅值与相位');
title('输入信号的频谱图');
grid on;
hold on;
%**********************************************************
%求输入信号的自相关函数
[Ry,lags]=xcorr(x1,'coeff');%求出自相关序列
lagst=lags/fs;
figure(2);%画出图形
plot((lagst),Ry,'k');
xlabel('tao');
ylabel('Ry');
title('输入信号的自相关函数');
grid on;
hold on;
%************************************%
%求输入信号功率谱
Sqy=fft(Ry,length(Ry));%对自相关函数进行傅里叶变换,得到均方根谱
f2=(0:length(Sqy)-1)'*fs/length(Sqy);
magSqy=abs(Sqy);
figure(3);
plot(f2,magSqy,'r');
xlabel('频率(Hz)');
ylabel('功率谱密度');
title('输入信号功率谱');
grid on;
hold on;
%**********************************************
%生成测试低通滤波器
fp=2500;fq=3500;Fs=44100;
rp=1;rs=35;
wp=2*pi*fp/Fs;
ws=2*pi*fq/Fs;
wap=tan(wp/2);
was=tan(ws/2);
[n,Wn]=buttord(fp/Fs,fq/Fs,rp,rs,'s');
[b,a]=butter(n,Wn);
[z,p,k]=buttap(n);
[bp,ap]=zp2tf(z,p,k);
[bs,as]=lp2lp(bp,ap,wap);
[bz,az]=bilinear(bs,as,0.5);
[h,w]=freqz(bz,az,512,Fs);
figure(4)
plot(w,abs(h)),axis([0 5000 0 1.1]);
grid on; %测试波器滤器
title('滤波器频率特性')
%******************************
%通过低通滤波器
[k,l] = butter(n,Wn);
Y=filter(k,l,x1); %信号通过低通滤波器
figure(5);
plot(t,Y);
xlabel('时间(t)');
ylabel('幅值')
title('滤波后时域谱(巴氏滤波器)');
grid on;
%************************************
%两重自相关提取微弱信号
z=xcorr(Y,'unbiased');
[x2,lagss]=xcorr(z,'unbiased');
tt=lagss/fs;
figure(6)
plot(tt,x2,'b');
grid on;
title('两重自相关波形')
x22=zeros(1,49);
for n=1:49
for j=0:20
x22(n)=x22(n)+Y(n+48*j);
end
end
x22=x22/14;
t2=0:48;
figure(10)
plot(t2/fs,x22,'r');
title('周期算法波形')
grid on;
x2=x22;%有三角函数自相关函数和原函数关系知
%***********************************************************
m2=mean(x2); %求均值
v2=var(x2); %求方差
w2=v2+m2^2; %求均方值
%******************************************
%进行FFT变换并做频谱图
y2=fft(x2,N);%进行fft变换
magy2=abs(y2);%求幅值
angley2=angle(y2);%求相位
fy2=(0:length(y2)-1)*fs/length(y2);%进行对应的频率转换
figure(7);
plot(fy2,magy2,'k',fy2,angley2,'r');%做频谱图
xlabel('频率(Hz)');
ylabel('幅值与相位');
title('输出信号的频谱图');
grid on;
hold on;
%**********************************************************
%求输出信号的自相关函数
[Ry2,lags2]=xcorr(x2,'coeff');%求出自相关序列
lagst2=lags2/fs;
figure(8);%画出图形
plot((lagst2),Ry2,'k');
xlabel('tao');
ylabel('Ry');
title('输出信号的自相关函数');
grid on;
hold on;
%************************************%
%求输出信号功率谱
fc2=fft(Ry2);
cm2=abs(fc2);
f2=(0:length(fc2)-1)'*fs/length(fc2);
figure(9)
plot(f2,cm2,'r');
xlabel('频率(Hz)');
ylabel('功率谱密度');
title('输出信号功率谱');
grid on;
hold on;
%**********************************************