数字信号处理实验二：频率采样型FIR滤波器结构特点与实现

一、实验目的

1、学习使用频率采样型结构实现FIR滤波器，初步熟悉FIR滤波器的线性相位特点。

2、直观体会频率采样型滤波器所具有的“滤波器组”特性，即在并联结构的每条支路上可以分别得到输入信号的各次谐波。

3、学习使用周期冲激串检测所实现滤波器的频域响应。

二、实验内容

        1.周期N=16

        2.由幅频特性可得，信号具有直流分量、基频、二次谐波和三次谐波，因此仅有第0、1、2、3、13、14、15这几条谱线非0，且幅度与构造一值。相频特性奇对称，符合正弦信号两根谱线相反的特点。

        3.H(e^jw)符合低通特性，其在 H(k)处取值为精确值，其他点处为经内插函数叠加得到。

        4.输出信号中已没有三次谐波分量, 且各谐波的相位也依次发生了变化。

        5.第1路波形是直流分量的输出，因H(0)=1，幅度没有改变。第2、3路分别是基波以及二次谐振输出的波形。第4路信号波形为零，因为H(3)=0，三次谐波被滤波器滤掉了。

        6.冲激串信号含有全部频率，故输出信号的时域特性和第三步中h(n)相同，幅频特性和第三步中的幅频特性在特定点的采样，近似为窗型。

三、实验代码

1.

clc;clear;

f=50;
N=16;
fs=N*f; 	%时域信号的采样频率
T=1/fs;		%采样间隔
L=N*2; 		%采样点数
t=0:T:(L-1)*T;		%构造输入信号

s0=0.5*cos(0);
s1=1*cos(2*pi*f*t+pi/2);
s2=0.5*cos(2*pi*2*f*t+pi);
s3=2*cos(2*pi*3*f*t-pi/2);
s=s0+s1+s2+s3;

i=0:1:L-1;
stem(i,s(i+1));
title('时域采样结果');
grid on;

 2.

clc;clear;

f=50;
N=16;
fs=N*f; 	%时域信号的采样频率
T=1/fs;		%采样间隔
L=N*2; 		%采样点数
t=N*T:T:(L-1)*T;

s0=0.5*cos(0);
s1=1*cos(2*pi*f*t+pi/2);
s2=0.5*cos(2*pi*2*f*t+pi);
s3=2*cos(2*pi*3*f*t-pi/2);
s=s0+s1+s2+s3;

a=fft(s);	%FFT
y=abs(a);   %幅度
z=angle(a); %相位
i=0:1:N-1;

subplot(2,1,1);
stem(i,y(i+1));
title('幅频特性');
subplot(2,1,2);
stem(i,z(i+1));
title('相频特性');

3.

clc;clear;

N=16;
H = zeros(1,N); 
H(1) = 1; 
H(2) = exp(-1i * pi * (N-1)/N);
H(3) = exp(-1i * 2 * pi * (N-1)/N);
H(15) = -exp(-1i * 14 * pi * (N-1)/N);
H(16) = -exp(-1i * 15 * pi * (N-1)/N);

% 抽头系数
h = ifft(H);    %ifft
subplot(2,2,1);
stem(0:N-1,real(h));
title('h(n)');

% 插值函数
k = 1:N;
Hz = zeros(1,N);    % H(Z)频率特性
a = 0;
for w = 0:0.005:2*pi
 a = a + 1;
 z = exp(1i*w);
 Hz(a) = 1/N*(1 - power(z,-N))*sum(H(k)./(1 - exp(1i*2*pi/N*(k-1)) / z ));%频率响应值
end

% 幅频特性
subplot(2,2,3);
plot(0:0.005/pi:2,abs(Hz))
title('H(z)幅频特性');

% 相频特性
subplot(2,2,4);
plot(0:0.005/pi:2,angle(Hz))
title('H(z)相频特性');

4.

clc;clear;
%输入信号
f=50;
N=16;
fs=N*f;
L=2*N;
T=1/fs;
t=0:T:(L-1)*T;
s0=0.5*cos(0);
s1=1*cos(2*pi*1*f*t+pi/2);
s2=0.5*cos(2*pi*2*f*t+pi);
s3=2*cos(2*pi*3*f*t-pi/2);
s=s0+s1+s2+s3;

x=s;
r=0.999;
y=CombFilter(x,N,r);    %对输入信号x滤波

x=y;
z=zeros(1,48);
for i=0:1:(N/2)
	Order=i;    %滤波器阶数
	H=[1 exp(-1i*pi*(N-1)/N) exp((-1i*2*pi*(N-1)/N)) 0 0 0 0 0 0 0 0 0 0 0 -exp((-1i*14*pi*(N-1)/N)) -exp((-1i*15*pi*(N-1)/N))];
	y=Resonator2(x,N,r,Order,H(i+1));
	z=z+y;
end
y=z/N;
for i=N:1:L-1
    y1(i-N+1)=y(i+1);  %取出N到L-1的元素
end

%时域
i=0:1:N-1;
subplot(2,2,1),stem(i,y1(i+1));
xlabel('n');
ylabel('振幅');
title('时域波形');
grid on;
 
%频谱
y=fft(y1);
y1=abs(y);
ang=angle(y);
i=0:1:N-1;

subplot(2,2,3),stem(i,y1(i+1));
title('幅频特性');
grid on;

subplot(2,2,4),stem(i,ang(i+1));
title('相频特性');
grid on;

5.

clc;clear;

f=50;
N=16;
fs=N*f; 	
T=1/fs;		
L=N*2; 		
t0=0:T:(L-1)*T;
s0=0.5*cos(2*pi*0*f*t0+0);
s1=1*cos(2*pi*f*t0+pi/2);
s2=0.5*cos(2*pi*2*f*t0+pi);
s3=2*cos(2*pi*3*f*t0-pi/2);
s=s0+s1+s2+s3;			%生成信号

%前4路谐振器的输出
r=0.999;
y1=CombFilter(s,N,r);
H=[1,exp(-1i*pi*(N-1)/N),exp(-1i*2*pi*(N-1)/N),0,0,0,0,0,0,0,0,0,0,0,exp(-1i*14*pi*(N-1)/N),exp(-1i*15*pi*(N-1)/N)];

y_h0=Resonator2(y1,N,r,0,H(1));
y_h0=y_h0/N;
y_h0=y_h0(N+1:2*N);

y_h1=Resonator2(y1,N,r,1,H(2));
y_h1=y_h1/N;
y_h1=y_h1(N+1:2*N);

y_h2=Resonator2(y1,N,r,2,H(3));
y_h2=y_h2/N;
y_h2=y_h2(N+1:2*N);

y_h3=Resonator2(y1,N,r,3,H(4));
y_h3=y_h3/N;
y_h3=y_h3(N+1:2*N);

%输出图像
k=0:N-1;
figure(1);

subplot(2,2,1);
stem(k,y_h0);
title('第一路谐振器输出');
subplot(2,2,2);
stem(k,y_h1);
title('第二路谐振器输出');
subplot(2,2,3);
stem(k,y_h2);
title('第三路谐振器输出');
subplot(2,2,4);
stem(k,y_h3);
title('第四路谐振器输出');


figure(2);

t0 = 0:(N-1); 
subplot(2,2,1);
stem(t0, s0(1:N));
title('输入信号第一路');

subplot(2,2,2);
stem(t0, s1(1:N));
title('输入信号第二路');

subplot(2,2,3);
stem(t0, s2(1:N));
title('输入信号第三路');

subplot(2,2,4);
stem(t0, s3(1:N));
title('输入信号第四路');

6.

clc;clear;

N=16;
x=zeros(1,2*N);
x(1)=1;
x(N+1)=1;	%冲激串
H=zeros(1,N);
h=zeros(1,N);
H=[1,exp(-1i*pi*(N-1)/N),exp((-1i*2*pi*(N-1)/N)),0,0,0,0,0,0,0,0,0,0,0,-exp((-1i*14*pi*(N-1)/N)),-exp((-1i*15*pi*(N-1)/N))];
r=0.999;
 
y=CombFilter(x,N,r);	%冲激串滤波后
y1=zeros;
for  k=0:N/2
    y1=y1+Resonator2(y,N,r,k,H(k+1));
end

subplot(2,2,1)
k = 0:1:2*N-1;
stem(k, x)
title('输入信号时域')
grid on

y0=y1/N;    %经谐振器后
k=N:2*N-1;
hold on;
subplot(2,2,2);
stem(k,y0(k+1));
title('输出信号时域');
grid on;

Y=zeros(1,N); 
for k=N:2*N-1
Y(k-N+1)=y0(k+1);	%输出信号第二个周期
end
F=fft(Y);

k=0:1:N-1;
hold on
subplot(2,2,3)
stem(k,abs(F(k+1)))
title('幅频特性');
grid on;
subplot(2,2,4)
stem(k,angle(F(k+1)))
title('相频特性');