多相滤波器的实现

变采样和多相滤波器的实现。本程序实现了一个变采样程序，中间设计滤波器设计和插值抽取。其中滤波器设计用的是窗函数法，根据要求设计窗函数，得到窗函数的长度。接着是插值，滤波，抽取，得到最后变采样之后的波形文件、另外对比了用直接卷积和多相分解卷积两种方法最后的结果。

 clc;
 close all;
 clear all;
 load（‘InputData.mat‘）;
 
 x = InputData;
 
 % let us see it;
 plot（x）;
 % let us save it
 fs_in =  2.285714e+03 ;
 WAVEFILE = ‘input.wav‘;
 wavwrite（xfs_inWAVEFILE）;
 fs_out= 8000;
 % let us listen it
 sound（xfs_in）;
 
 % output and input ratio is
 disp（[‘Sameprate Convert Ratio is ‘num2str（fs_out/fs_in）]）;
 
 %3.5 =7/2
 disp（‘find a couple of relatively prime numbers 7/2 =3.5‘）;
 U=7;D=2;
 
 fprintf（‘Upsample ratio is %d\n‘U）;
 fprintf（‘Downsample ratio is %d\n‘D）;
 
 %Part A
 % transition band = 0.157080;过渡带 （Transition-band = pi/20）
 % pass band attenuation <= 7.000000 dB 
 % stop band attenuation >= 40.000000 dB 
 % comparing Anti-alias filter and Anti mirror filter
 % Wc = min（pi/7pi/2）Wc = pi/7
 % 用海明窗函数设计低通滤波器
 wc = pi/7;
 ws=wc+pi/40;
 wp=wc-pi/40;
 
 %Part B
 deltaw= ws - wp;        % 过渡带宽Δω的计算
 N = ceil（6.6*pi/ deltaw） + 1;    % 按海明窗计算所需的滤波器阶数N0
 wdham = （hamming（N+1））‘;        % 海明窗计算
 Wn=wc/（pi）      ;             % 计算截止频率这里要用归一化的，否则fir1函数不认
 b=fir1（NWnwdham）;
 [Hw]=freqz（b11024）;
 db=20*log10（abs（H））;
 % 画频响曲线
 figure
 fs_scaled =1;
 plot（w*fs_scaleddb）;title（‘幅度响应（单位： dB）‘）;grid
 figure
 fs_scaled =1;
 plot（w*fs_scaleddb）;title（‘幅度响应（单位： dB）‘）;grid
 axis（[0 pi -45 5]）; xlabel（‘归一化频率‘）; ylabel（‘分贝‘）;
 set（gca‘XTickMode‘‘manual‘‘XTick‘[0pi]）;
 set（gca‘YTickMode‘‘manual‘‘YTick‘[-455]）;

 %系统能承受1024点fft，可以认为把滤波器设置为这么长的长度
 delta_w= pi/1024;                            %搜索;滤波器是1024个点，对应着pid的位置，对应着fs/2的位置
 Ap=-（min（db（1:1:wp/delta_w+1））） ;  %实际通带纹波
 As=-round（max（db（ws/delta_w+1:1:513）））  ;  %实际阻带纹波
 
 fprintf（‘Actuall pass attenuation %f dB required <=7dB\n‘Ap）;
 fprintf（‘Actuall stop attenuation %f  dB required >=40dB\n‘As）;
 
%   Part C
% 请清参考附图

%  Part D 
%  首先进行7倍的插值
 U=7;
  
 Norignal=length（x）;
 %  Norignal=1024;
 x=x（1:Norignal）; %截取前N个点计算不全算的话可以取前1024个点
 x1=zeros（1U*length（x））;   
 x1（1:U:end）=x;   %对信号进行L倍直接插值
 
N = length（x1）;
%多相分解 134 阶滤波器 134 = 67 * 2
for r=1:67  %每相67个点
    for k=1:2  %2倍抽取
        hh（kr）=b（（r-1）*2+k）;
    end
end
%各路信号时延
for k=1:2
    for m=1:N-2
        xk（km）=x1（m+k-1）;  %信号进入各相经过不同的延时
    end
end
%对时延后的信号2倍抽取
for k=1:2
    for m=1:2:N-2
        xd（k（m-1）/2+1）=xk（km）;
    end
end
%多相滤波
n2=N/2-1;  %抽取后信号的点数
for k=1:2
    for m=1:n2
        xdk（m）=xd（km）;
    end
    for m=1:67
        hhk（m）=hh（km）;
    end
    y=conv（hhkxdk）;n3=n2+67-1;  %卷积之后信号的点数=N/2-1+14-1=212
    for m=1:n3
        y1（km）=y（m）;
    end
end
yt=y1（1:）+y1（2:）  ;%合并两相

% 直接做卷积的结果
% 7倍抽取之后直接送入LPF
x2  = conv（bx1）;
% 再对x2做两倍的抽取得到x3
for r=1:length（x2）/2  
    x3（r） = x2（（r-1）*2+1）;
end
% x3就是直接滤波结果，把两者作比较：
figure
plot（[100:200]yt（100:200）‘b‘）;
hold on
plot（[100:200]x3（100:200）‘r‘）;
title（‘Output comparasion by 2 ways‘）;
legend（‘polyphase‘‘conventional‘）;
% output the last wave files at 8k samperate
wavwrite（3.5*x3fs_out‘output1.wav‘）;
w

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件     906966  2014-12-13 11:53  PolyphaseFilter\block diagram for a polyphase decomposition.bmp

     文件      14335  2014-12-13 15:40  PolyphaseFilter\Design.JPG

     文件       2977  2014-12-11 09:40  PolyphaseFilter\Email.txt

     文件       3948  2014-12-13 15:17  PolyphaseFilter\FilterDesign.asv

     文件       3859  2014-12-13 15:20  PolyphaseFilter\FilterDesign.m

     文件      13760  2014-12-13 15:07  PolyphaseFilter\input.wav

     文件      52974  2014-12-11 09:38  PolyphaseFilter\InputData.mat

     文件      48182  2014-12-13 14:33  PolyphaseFilter\output1.wav

     文件      48180  2014-12-13 14:33  PolyphaseFilter\output2.wav

     文件      70386  2014-12-13 10:19  PolyphaseFilter\六种窗函数特性表.jpg

     文件        139  2014-12-13 11:52  PolyphaseFilter\参考链接.txt

     目录          0  2016-09-20 15:52  PolyphaseFilter

----------- ---------  ---------- -----  ----

              1165706                    12