插值滤波相关文档(含matlab代码)

上篇文字里面，我们说到了FIR滤波器，接下来我们在此处讨论一下插值滤波。

这个话题属于多速率信号处理的范畴，经典的参考读物有杨小牛老师的《软件无线电原理与应用》以及Fredric J. Harris 所著的《通信系统中的多采样率信号处理》，这两位都是通信行业中的大腕，前者曾是中电科36所的所长，这个所在嘉兴，是国内电子对抗领域的翘楚。后面一位更加厉害，高通公司用他的名字在圣迭戈州立大学设立了一个奖学金。

数字上、下变频是软件无线电领域里面非常重要的一个分支，以数字上变频来说，关键技术就是插值滤波和DDS技术。插值滤波技术，也叫做上采样，或采样率提升，为什么这个技术很重要呢，请自行搜索“软件无线电插值滤波器”。

根据小牛老师的书，信号N倍插值之后的频谱会出现N个镜像。于是呢，需要用滤波器来把镜像频率除去，插值滤波器常用的有2个系列，FIR和CIC，其中FIR主要用于小倍率的插值，而CIC则用于大倍率的插值，比如说，参照以下ADI 的经典芯片AD9857，里面用FIR做了4倍插值，用CIC做2-63倍插值，这是因为FIR需要进行大量乘加，但是滤波性能好，而CIC不用乘法，但是只能抑制镜像频率附近的信号频带，并且CIC的通带内衰减较快，所以CIC适合用在高采样率的一端，关于CIC先不多说，此处专注FIR系列。

贴一张小牛老师书里面的插图

多速率FIR滤波器有两种说法，半带滤波器（halfband filter）和多相滤波器，半带其实也是多相的一个特例，区别在于，如果进行2的N次幂的插值，我们通常用半带，因为乘法数量可以降为原来的1/4，而对于非基2的插值，我们就得用多相结构了。无论是半带还是多相，其实都是一种对原型滤波器的优化。噢，原型滤波器，对了，为了设计最终使用的滤波器，我们首先要有一个原型滤波器，然后根据这个原型滤波器再设计半带或是多相结构的滤波器，详情请自行拜读小牛爷爷的书吧，祖师爷的书最好还是买一本算是交了保护费嘛。

此处的仿真过程大致如下，首先，生成一个多音正弦，然后进行补零的插值，然后把补零插值的信号塞到一个低通滤波器里面，这样就有三个版本的信号，原始信号，补零插值信号，以及插值滤波之后的信号。注意此仿真的计算过程是概念上的理论模型，实际干活中不会这么做的，因为进入滤波器的数据里面包含的大量的0数据，把这种数据拿去做乘法是个非常亏本的事情，所以才会有半带啊，

多相啊这种优化，嗯，不再赘述了，先看图吧，另外，吐槽一下，EDN网站啊，你是我们码农电工的乐园，尽管杂志是免费送的，网站是免费看到，博客是免费写的，您就不能支持一下清楚一点的图片么，您瞧我这几张时域频域信号图片寒碜的，我自己都觉得难受，实在不成您办个募捐，俺们捐点款给您买带宽买磁盘阵列也成啊。

三个版本信号的时域图

滤波器频响

原始信号频谱补零插值信号频谱

插值滤波信号频谱OK，上代码，请自行修改参数绘图，祝玩得开心

%///////////////////////////////////////////////////////////

% FILE: test1.m

% Simulation of FIR interpolation

%///////////////////////////////////////////////////////////

close all;

clear;

clc;

% specify input multi-tone sine component frequency

sin_freq= [1:3]*20E3;

data_len= 2048;% signal data length

fs= 600E3 ;% sample rate

quant_bits= 12;% signal quant bits

kaiser_beta = 8;% beta of kaiser win

n_itp= 3;% num of interpolation

x_max_1x= 30;% 1x signal time plot x axis max value

n_coef= 32;% number of coefficents

f= [0, 0.2, 0.2, 1];

m= [1, 1, 0, 0];

coeff = fir2(n_coef-1, f, m);

% draw the filter reponse curve

freqz(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 signal

signal_1x= gen_quant_multi_sin(fs, sin_freq, data_len, quant_bits); % create zero filling interpolation signal

signal_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);

filter_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 ', 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]);% 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 range

title('Zero Interpolation, Normalized Spectrum', 'fontsize', 14);

kaiser_win_spectrum_plot(fs_itp, filter_out, kaiser_beta);

ylim([-160,10]);% set y-axis range

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