插值濾波相關(guān)文檔(含matlab代碼)

1、上篇文字里面，我們說(shuō)到了FIR濾波器，接下來(lái)我們?cè)诖颂幱懻撘幌虏逯禐V波。這個(gè)話(huà)題屬于多速率信號(hào)處理的范疇，經(jīng)典的參考讀物有楊小牛老師的軟件無(wú)線(xiàn)電原理與應(yīng)用以及 Fredric J. Harris 所著的 通信系統(tǒng)中的多采樣率信號(hào)處理，這兩位都是通信行業(yè)中的大腕，前者曾是中電科36所的所長(zhǎng)，這個(gè)所在嘉興，是國(guó)內(nèi)電子對(duì)抗領(lǐng)域的翹楚。后面一位更加厲害，高通公司用他的名字在圣迭戈州立大學(xué)設(shè)立了一個(gè)獎(jiǎng)學(xué)金。數(shù)字上、下變頻是軟件無(wú)線(xiàn)電領(lǐng)域里面非常重要的一個(gè)分支，以數(shù)字上變頻來(lái)說(shuō)，關(guān)鍵技術(shù)就是插值濾波和DDS技術(shù)。插值濾波技術(shù)，也叫做上采樣，或采樣率提升，為什么這個(gè)技術(shù)很重要呢，請(qǐng)自行搜索“軟件無(wú)線(xiàn)電 插值

2、濾波器”。根據(jù)小牛老師的書(shū)，信號(hào)N倍插值之后的頻譜會(huì)出現(xiàn)N個(gè)鏡像。于是呢，需要用濾波器來(lái)把鏡像頻率除去，插值濾波器常用的有2個(gè)系列，F(xiàn)IR和CIC，其中FIR主要用于小倍率的插值，而CIC則用于大倍率的插值，比如說(shuō)，參照以下ADI的經(jīng)典芯片 AD9857，里面用FIR做了4倍插值，用CIC做2-63倍插值，這是因?yàn)镕IR需要進(jìn)行大量乘加，但是濾波性能好，而CIC不用乘法，但是只能抑制鏡像頻率附近的信號(hào)頻帶，并且CIC的通帶內(nèi)衰減較快，所以CIC適合用在高采樣率的一端，關(guān)于CIC先不多說(shuō)，此處專(zhuān)注FIR系列。貼一張小牛老師書(shū)里面的插圖多速率FIR濾波器有兩種說(shuō)法，半帶濾波器（halfband f

3、ilter）和多相濾波器，半帶其實(shí)也是多相的一個(gè)特例，區(qū)別在于，如果進(jìn)行2的N次冪的插值，我們通常用半帶，因?yàn)槌朔〝?shù)量可以降為原來(lái)的1/4，而對(duì)于非基2的插值，我們就得用多相結(jié)構(gòu)了。無(wú)論是半帶還是多相，其實(shí)都是一種對(duì)原型濾波器的優(yōu)化。噢，原型濾波器，對(duì)了，為了設(shè)計(jì)最終使用的濾波器，我們首先要有一個(gè)原型濾波器，然后根據(jù)這個(gè)原型濾波器再設(shè)計(jì)半帶或是多相結(jié)構(gòu)的濾波器，詳情請(qǐng)自行拜讀小牛爺爺?shù)臅?shū)吧，祖師爺?shù)臅?shū)最好還是買(mǎi)一本算是交了保護(hù)費(fèi)嘛。此處的仿真過(guò)程大致如下，首先，生成一個(gè)多音正弦，然后進(jìn)行補(bǔ)零的插值，然后把補(bǔ)零插值的信號(hào)塞到一個(gè)低通濾波器里面，這樣就有三個(gè)版本的信號(hào)，原始信號(hào)，補(bǔ)零插值信號(hào)，以及

4、插值濾波之后的信號(hào)。注意此仿真的計(jì)算過(guò)程是概念上的理論模型，實(shí)際干活中不會(huì)這么做的，因?yàn)檫M(jìn)入濾波器的數(shù)據(jù)里面包含的大量的0數(shù)據(jù)，把這種數(shù)據(jù)拿去做乘法是個(gè)非常虧本的事情，所以才會(huì)有半帶啊，多相啊這種優(yōu)化，嗯，不再贅述了，先看圖吧，另外，吐槽一下，EDN網(wǎng)站啊，你是我們碼農(nóng)電工的樂(lè)園，盡管雜志是免費(fèi)送的，網(wǎng)站是免費(fèi)看到，博客是免費(fèi)寫(xiě)的，您就不能支持一下清楚一點(diǎn)的圖片么，您瞧我這幾張時(shí)域頻域信號(hào)圖片寒磣的，我自己都覺(jué)得難受，實(shí)在不成您辦個(gè)募捐，俺們捐點(diǎn)款給您買(mǎi)帶寬買(mǎi)磁盤(pán)陣列也成啊。三個(gè)版本信號(hào)的時(shí)域圖濾波器頻響原始信號(hào)頻譜補(bǔ)零插值信號(hào)頻譜插值濾波信號(hào)頻譜OK，上代碼，請(qǐng)自行修改參數(shù)繪圖，祝玩得開(kāi)心%

5、/% FILE: test1.m% Simulation  of FIR interpolation%/close all;clear    ;clc      ;% specify input multi-tone sine component frequencysin_freq    = 1:3*20E3;data_len    = 2048  ;% signal data lengthfs   

6、60;      = 600E3 ;% sample ratequant_bits  = 12    ;% signal quant bitskaiser_beta = 8     ;% beta of kaiser winn_itp       = 3     ;% num of interpolationx_max_1x    = 3

7、0    ;% 1x signal time plot  x axis max valuen_coef      = 32    ;% number of coefficentsf     = 0, 0.2, 0.2, 1    ;m     = 1, 1  , 0  , 0    ;coeff = fir2(n_co

8、ef-1, f, m);% draw the filter reponse curvefreqz(coeff);x_max_itp   = x_max_1x * n_itp ;fs_itp      = fs * n_itp       ;% itp for abv of 'interpolation'data_len_itp= data_len * n_itp ;% generate original signalsignal_1x &#

9、160; = gen_quant_multi_sin(fs, sin_freq, data_len, quant_bits);% create zero filling interpolation signalsignal_itp  = zeros(data_len_itp,1);% write the original signal value into itp signal signal_itp(1:n_itp:data_len_itp-n_itp+1) = signal_1x;data_conv  = conv(coeff, signal_itp);filt

10、er_out = data_conv(1:data_len_itp);figure;subplot(3,1,1);stem(signal_1x , 'MarkerSize', 1.5); xlim(1 x_max_1x);tt_str = 'Signal 1x' title(tt_str , 'fontsize', 14);subplot(3,1,2);stem(signal_itp, 'MarkerSize', 1.5);xlim(1 x_max_itp);tt_str = strcat('Signal Itp

11、 ', num2str(n_itp), 'x'); title(tt_str, 'fontsize', 14);subplot(3,1,3);stem(filter_out, 'MarkerSize', 1.5);xlim(1 x_max_itp);tt_str = 'Filter Out' title(tt_str, 'fontsize', 14);kaiser_win_spectrum_plot(fs, signal_1x, kaiser_beta);ylim(-160,10)

12、60;    ;   % set y-axis range        title('signal 1x, Normalized Spectrum', 'fontsize', 14);kaiser_win_spectrum_plot(fs_itp, signal_itp, kaiser_beta);ylim(-160,10)     ;   % set y-axis ran

13、ge        title('Zero Interpolation, Normalized Spectrum', 'fontsize', 14);kaiser_win_spectrum_plot(fs_itp, filter_out, kaiser_beta);ylim(-160,10)     ;   % set y-axis range       

14、; title('Filter Out, Normalized Spectrum','fontsize', 14);% % CIC compensating filter design using frequency sampling method clear all close all % CIC filter parameters % R = 4; % Decimation factor M = 2; % Differential Delay N = 4; % Number of Stages B = 17; % Number of bits to

15、 represent fixed point filter coefficients Fs = 80e6; % (High) Sampling frequency in Hz (before decimation) Fc = 4e6; % Passband edge in Hz % fir2.m parameters % L = 30; % Order of filter taps; must be an even number Fo = R*Fc/Fs; % Normalized Cutoff freq; 0<Fo<=0.5/M; % Fo should be less than

16、 1/(4M) for good performance % Fo = 0.5/M; % use Fo=0.5 if we don't care responses outside passband % CIC Compensator Design using fir2.m % p = 2e3; % Granulatiry s = 0.25/p; % Stepsize fp = 0:s:Fo; % Passband frequency samples fs = (Fo+s):s:0.5; % Stopband frequency samples f = fp fs*2; % Noram

17、lized frequency samples; 0<=f<=1; Mp = ones(1,length(fp); % Passband response; Mp(1)=1 Mp(2:end) = abs( M*R*sin(pi*fp(2:end)/R)./sin(pi*M*fp(2:end).N; % Inverse sinc Mf = Mp zeros(1,length(fs); f(end) = 1; h = fir2(L,f,Mf); % Filter length L+1 h = h/norm(h); % Floating point coefficients hz =

18、floor(h*power(2,B); % Quantization of filter coefficients % Full resolution CIC filter response % hrec = ones(1,R*M); tmph = hrec; for k=1:N-1 tmph = conv(hrec, tmph); end; hcic = tmph; hcic=hcic/norm(hcic); % Total Response % hzp = upsample(hz,R); hp = upsample(h, R); ht = conv(hcic, hp); % Concate

19、nation of CIC and fir2 FIR at high freqency hzt = conv(hcic, hzp); % CIC + Fixed point fir2 at high frequency Hcic, wt = freqz(hcic, 1, 4096, Fs); % CIC Freq. Response Hciccomp, wt = freqz(hp, 1, 4096, Fs); % CIC Comp. response using fir2 Ht, wt = freqz(ht, 1, 4096, Fs); % Total response for CIC + f

20、loating point fir2 Hzt, wt = freqz(hzt, 1, 4096, Fs); % Total response for CIC + fixed point fir2 Mcic = 20*log10(abs(Hcic)/max(abs(Hcic); % CIC Freq. Response Mciccomp = 20*log10(abs(Hciccomp)/max(abs(Hciccomp); % CIC Comp. response using fir2 Mt = 20*log10(abs(Ht)/max(abs(Ht); % Total response for CIC + floating point fir2 Mzt = 20*log10(abs(Hzt)/max(abs(Hzt); % Total response for CIC