qmf.asv

这是我的论文要设计程序,通过设计QMF(CQF,OQF)滤波器,研究针对不同的滤波器阶数

clear
format long g;
% n=input('Please Input the Number of the Filter:');%
fc=input('Please Input thecut-off frequency in 4950--5050（For PR）:');
for n = 6 : 100
fs=10000;%采样频率
%fc=4965;%为使其滤波器响应接近最接近１
wn=fc/fs;
b1=fir1(n,wn,'stop',rectwin(n+1),'noscale');   %
%b1=[b1,zeros(1,250-length(b1))];%补零(圆周卷积,且增大分辨率便于观察)
%%QMFl2 ：“l”代表滤波，“２“代表第二通道；QMFｃ２：“ｃ”代表重建，“２”代表第二通道
for k=1:n+1   QMFl2(k)=((-1)^(k-1))*b1(k);%(k-1)是因为这里ｋ从１开始不是从０开始；
              QMFc1(k)=b1(k); 
              QMFc2(k)=(-1)^(k)*b1(k); 
              CQFl2(k)=((-1)^(k-n-1))*b1(n-k+2); 
              CQFc2(k)=QMFc2(k);
              OQFl2(k)=(-1)^(n+1-k)*b1(n+2-k);
              OQFc1(k)=b1(n+2-k);
              OQFc2(k)=(-1)^(k-1)*b1(k);
end
for k=1:n+1 
    CQFc1(k)=(-1)^(k-1)*CQFl2(k);        %CQF第一通道重建系数
end


 
[H1z,w]=freqz(b1,1);      h1=abs(H1z);   g1=20*log10(h1);
[H2z,w]=freqz(QMFl2,1);   h2=abs(H2z);   g2=20*log10(h2);
[H3z,w]=freqz(CQFl2,1);   h3=abs(H3z);   g3=20*log10(h3);
[H4z,w]=freqz(OQFl2,1);   h4=abs(H3z);   g4=20*log10(h4);

%%%%%画滤波器的频谱图，相位图
figure(1)
subplot(131);
plot(w/pi,g1,'g',w/pi,g2,'b') %QMF幅频图
title('QMF幅频图');

subplot(132)
plot(w/pi,g1,'g',w/pi,g3,'b') %CQF 幅频图
title('CQF幅频图');
subplot(133);
plot(w/pi,g1,'g',w/pi,g4,'b');
title('OQF幅频图');

%%%%%输入信号滤波
t=(100:1000)/fs;
inpfilenm = 's2ofwb'; 
%xt=sin(2*pi*100*t)+sin(2*pi*400*t)+sin(2*pi*700*t);
[xt, fsingnal] =wavread(inpfilenm);     
Y=fft(xt,512);
Pyy=Y.*conj(Y)/512;
Pym=max(Pyy);
f=10000*(0:256)/512;
figure(2);
plot(f,Pyy(1:257)/Pym)
title('输入信号频谱图'); axis([0 1000 0 1]); axis manual; 


shuchu=filter(b1,1,xt); 

sshuchu=filter(QMFl2,1,xt); %QMF滤波

sshuchuc=filter(CQFl2,1,xt);%  CQF滤波
                            
sshuchuo=filter(OQFl2,1,xt);%OQF滤波器
                             
%%%%%%%%%%%%%%抽取

shuchu=dyaddown(shuchu); 
    
sshuchu=dyaddown(sshuchu); %QMF
    
sshuchuc=dyaddown(sshuchuc); %QMF抽取
                             
sshuchuo=dyaddown(sshuchuo);
                               
 %%%%%%%%%插值
   shuchu=dyadup(shuchu) ;
   
   sshuchu=dyadup(sshuchu);%CQF
  
   sshuchuc=dyadup(sshuchuc);%CQF
                               
   sshuchuo=dyadup(sshuchuo);
                               
shuchu1=filter(QMFc1,1,shuchu) ; 

sshuchu1=filter(QMFc2,1,sshuchu);%QMF重建

shuchuc1=filter(CQFc1,1,shuchu);   

sshuchuc1=filter(CQFc2,1,sshuchuc);%QMF重建

shuchuo1=filter(OQFc1,1,shuchu); 

sshuchuo1=filter(OQFc2,1,sshuchuo);

                     
shuchuhe=(shuchu1+sshuchu1)*2;%QMF输出

shuchuhec=(shuchuc1+sshuchuc1)*2;%CQF输出

shuchuheo=(shuchuo1+sshuchuo1)*2;

%figure(3);  plot(t-n/fs,shuchuhe,'g--',t,xt,'b'); axis([0.01 0.1-n/fs -5 5]);  legend('重构信号','原始信号');title('QMF重建与原始信号');     
%figure(4);  plot(t-n/fs,shuchuhec,'g--',t,xt,'b'); axis([0.01 0.1-n/fs -5 5]);  legend('重构信号','原始信号');title('CQF重建与原始信号');       
%figure(5)  ;plot(t-n/fs,shuchuheo,'g--',t,xt,'b');axis([0.01 0.1-n/fs -5 5]); legend('重构信号','原始信号'); title('OQF重建与原始信号');



figure(6);subplot(311);
num=8192;                       
f=fs*(0:num/2-1)/num;                  
Y=fft(shuchuhe,num);                           
P=abs(Y);                                     
P1=max(P);                                     
plot(f,P(1:num/2)/P1);           
title('QMF重建幅频响应');                       
axis([0 1000 0 1]); axis manual;
subplot(312);
Y=fft(shuchuhec,num);                           
P=abs(Y);                                     
P1=max(P);                                     
plot(f,P(1:num/2)/P1);           
title('CQF重建幅频响应');                       
axis([0 1000 0 1]); axis manual;                                     
subplot(313);
Y=fft(shuchuheo,num);                           
P=abs(Y);                                     
P1=max(P);                                     
plot(f,P(1:num/2)/P1);           
title('OQF重建幅频响应');                      
axis([0 1000 0 1]); axis manual;