mallatpu.m

用MALLAT算法对信号分解,并给出各分解信号的频谱

clear;
%N=512;
%x=rand(1,N);
x=wavread('duanyu.wav');
%N=length(x);
N=25000;
x=x(1:N);
subplot(2,1,1);
plot(x);
title('原始信号图');
p=fft(x);
Ay=2*sqrt(p.*conj(p))/25000;
subplot(2,1,2);
ff=0:(25000/2-1);
plot(ff,Ay(1:25000/2))
title('原始信号频谱图');
axis([0,25000/2,0,0.01])
x=x';
[ld,hd]=wfilters('db4','d');
%%%%%%%%%%%%%%% di yi ci fen jie %%%%%%%%%%%%%%%%%%%
ld1=[ld,zeros(1,N-length(ld))];
hd1=[hd,zeros(1,N-length(hd))];
a1=ifft(fft(x).*fft(ld1));
d1=ifft(fft(x).*fft(hd1));
a1=dyaddown(a1);
d1=dyaddown(d1);
%%%%%%%%%%%%%%%%%%% di er ci fen jie %%%%%%%%%%%%%%%%%%%%%
ld2=[ld,zeros(1,N/2-length(ld))];
hd2=[hd,zeros(1,N/2-length(hd))];
a2=ifft(fft(a1).*fft(ld2));
d2=ifft(fft(a1).*fft(hd2));
a2=dyaddown(a2);
d2=dyaddown(d2);
%%%%%%%%%%%%%%%%%%%% di san ci fen jie %%%%%%%%%%%%%%
ld3=[ld,zeros(1,N/4-length(ld))];
hd3=[hd,zeros(1,N/4-length(hd))];
a3=ifft(fft(a2).*fft(ld3));
d3=ifft(fft(a2).*fft(hd3));
a3=dyaddown(a3);
d3=dyaddown(d3);
%%%%%%%%%%%%%%%%%%%%%%%%% di si ci fen jie %%%%%%%%%%%%%%%%%%%
ld4=[ld,zeros(1,N/8-length(ld))];
hd4=[hd,zeros(1,N/8-length(hd))];
a4=ifft(fft(a3).*fft(ld4));
d4=ifft(fft(a3).*fft(hd4));
a4=dyaddown(a4);
d4=dyaddown(d4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%画出分解后信号频谱图%%%%%%%%%%%%%%%
figure(2)
subplot(4,1,1);
ap1=fft(a4);
ap1=2*sqrt(a4.*conj(a4))/1562;
ff=0:(1562/2-1);
plot(ff,ap1(1:1562/2));
subplot(4,1,2);
dp4=fft(d4);
dp4=2*sqrt(d4.*conj(d4))/1562;
ff=0:(1562/2-1);
plot(ff,dp4(1:1562/2));
subplot(4,1,3);
dp3=fft(d3);
dp3=2*sqrt(d3.*conj(d3))/3125;
ff=0:(3125/2-1);
plot(ff,dp3(1:3125/2));
subplot(4,1,4);
dp2=fft(d2);
dp2=2*sqrt(d2.*conj(d2))/6250;
ff=0:(6250/2-1);
plot(ff,dp2(1:6250/2));
figure(3)
dp1=fft(d1);
dp1=2*sqrt(d1.*conj(d1))/12500;
ff=0:(12500/2-1);
plot(ff,dp1(1:12500/2));