信号分析与处理实验报告

华北电力大学

实验报告

|

|

实验名称FFT的软件实现实验（Matlab）

IIR数字滤波器的设计

课程名称信号分析与处理

|

|

专业班级：电气化1308 学生姓名：袁拉麻加

学号： 201301000827 成绩：

指导教师：杨光实验日期： 2015-12-17

快速傅里叶变换实验

一、实验目的及要求

通过编写程序，深入理解快速傅里叶变换算法（FFT）的含义，完成FFT和IFFT算法的软件实现。

二、实验内容

利用时间抽取算法，编写基2点的快速傅立叶变换（FFT）程序；并在FFT

程序基础上编写快速傅里叶反变换（IFFT）的程序。

三：实验要求

1、FFT和IFFT子程序相对独立、具有一般性，并加详细注释；

2、验证例6-4，并能得到正确结果。

3、理解应用离散傅里叶变换（DFT）分析连续时间信号频谱的数学物理基础。

四、实验原理：

a.算法原理

1、程序输入序列的元素数目必须为2的整数次幂，即N=2M,整个运算需要M

级蝶形运算；

2、输入序列应该按二进制的码位倒置排列，输出序列按自然序列排列；

3、每个蝶形运算的输出数据军官占用其他输入数据的存储单元，实现“即位运算”；

4、每一级包括N/2个基本蝶形运算，共有M*N/2个基本蝶形运算；

5、第L级中有N/2L个群，群与群的间隔为2L。

6、处于同一级的各个群的系数W分布相同，第L级的群中有2L-1个系数；

7、处于第L级的群的系数是(p=1,2,3，…….,2L-1)而对于第L级的

蝶形

运算，两个输入数据的间隔为2L-1。

b.码位倒置程序流程图

c. 蝶形运算程序流程图

五、程序代码与实验结果

a.FFT

程序：

%%

clear all;close all;clc;

%输入数据%

A=input('输入x(n)序列','s');

A=str2num(A);

% A=[1,2,-1,4]; %测试数据%

%%

%校验序列,%

n=length(A);

m=log2(n);

if (fix(m)~=m)

disp('输入序列长度错误，请重新输入!');

A=input('输入x(n)序列','s');

A=str2num(A);

else

disp('输入正确，请运行下一步')

end

%%

%码位倒置%

for k=0:n-1

for j=1:m %取M位的二进制数%

x1(j)=bitget(k,j); %倒取出二进制数%

end

x1=num2str(x1); %将数字序列转化为字符串%

y(k+1)=bin2dec(x1); %二进制序列转化为十进制数%

clear x1

end

for k=1:n

B(k)=A(y(k)+1); %时间抽取序列%

end

clear A

%%

%计算%

for L=1:m %分解为M级进行运算%

LE=2^L; %第L级群间隔为2^L%

LE1=2^(L-1); %第L级中共有2^(L-1)个Wn乘数，进行运算蝶运算的两数序号相隔LE1%

W=1;

W1=exp(-1i*pi/LE1);

for R=1:LE1 %针对第R个Wn系数进行一轮蝶运算，共进行LE1次%

for P=R:LE:n %每个蝶的大小为LE% Q=P+LE1;

T=B(Q)*W;

B(Q)=B(P)-T;

B(P)=B(P)+T;

end

W=W*W1;

end

end

B %输出X(k)%

%%

验证结果：

例6-4

b.IFFT

程序：

%%

clear all;close all;clc;

%输入数据%

A=input('输入X(k)序列','s');

A=str2num(A);

% A=[6,2+2i,-6,2-2i]; %测试数据%

%%

%校验序列,%

n=length(A);

m=log2(n);

if (fix(m)~=m)

disp('输入序列长度错误，请重新输入!');

A=input('输入x(n)序列','s');

A=str2num(A);

else

disp('输入正确，请运行下一步')

end

%%

%码位倒置%

for k=0:n-1

for j=1:m %取M位的二进制数%

x1(j)=bitget(k,j); %倒取出二进制数%

end

x1=num2str(x1); %将数字序列转化为字符串%

y(k+1)=bin2dec(x1); %二进制序列转化为十进制数%

clear x1

end

for k=1:n

B(k)=A(y(k)+1); %时间抽取序列%

end

clear A

%%

%计算%

for L=1:m %分解为M级进行运算%

LE=2^L; %第L级群间隔为2^L%

LE1=2^(L-1); %第L级中共有2^(L-1)个Wn乘数，进行运算蝶运算的两数序号相隔LE1%

W=1;

W1=exp(-1i*pi/LE1);

for R=1:LE1 %针对第R个Wn系数进行一轮蝶运算，共进行LE1次%

for P=R:LE:n %每个蝶的大小为LE%

Q=P+LE1;

T=B(Q)*W;

B(Q)=B(P)-T;

B(P)=B(P)+T;

end

W=W*W1;

end

end

B=conj(B); %取共轭%

B=B/n %输出x(n)%

验证结果：

六、实验心得与结论

本次实验借助于Matlab软件，我避开了用C平台进行复杂的复数运算，在一定程度上简化了程序，并添加了简单的检错代码，码位倒置我通过查阅资料，使用了一些函数，涉及到十-二进制转换，数字-文本转换，二-文本转换，相对较复杂，蝶运算我参考了书上了流程图，做些许改动就能直接实现。

通过本次实验，我对FFT的过程更加熟悉了，尤其是W N系数的特性，可以通过累乘来实现W系数的更新。

IIR数字滤波器的设计实验

一、实验目的

1.掌握实现模拟信号数字化的编程方法。

2.掌握实现时间抽取快速傅里叶变换（FFT）的编程方法。

3.掌握设计IIR数字滤波器的双线性变换法。

4.加深对信号分析与处理的理解。

二、实验内容

1.生成三个不同频率的正弦信号

1()

x t、2()

x t、3()

x t，将它们叠加成一个信号()

x t。用满足采样定理的采样频率对模拟信号()

x t进行采样，得到离散时间序列()

x n。

2.编制子程序实现FFT运算，并对()

x n进行FFT运算，观察它的频谱。3.设计满足已知的设计指标的IIR数字滤波器，并用它对()

x n进行滤波处理。对滤波后得到的波形再做FFT运算，观察滤波后的信号频谱。

三、实验结果

1.生成的信号

图1 原始信号及合成信号的时域波形

图2 原始信号的幅度谱

2.滤波器的性能检验

图3 模拟滤波器的频率特性

图4 数字滤波器的频率特性

3.滤波后的信号

图5 滤波前后信号的时域波形对比

图6 滤波前后信号的幅度谱对比

附录

附录一：FFT子程序

function [ A ] = fftlzr( A )

N=length(A);M=log2(N);

j=complex(0,1);

if fix(M)~=M

A=[A(:);zeros(2^(fix(M)+1)-N,1)]; end

N=length(A);M=log2(N);

I=0;J=0;L=1;

while I~=N-1

if I

t=A(J+1);

A(J+1)=A(I+1);

A(I+1)=t;

end

K=N/2;

while K<=J

J=J-K;

K=K/2;

end

J=J+K;

I=I+1;

end

while L<=M

LE=2^L; R=1;

LE1=LE/2;

W=1;W1=exp(-j*pi/LE1);

while R<=LE1

P=R;

while P<=N

Q=P+LE1;T=A(Q)*W; A(Q)=A(P)-T;

A(P)=A(P)+T;

P=P+LE;

end

W=W*W1;

R=R+1;

end

L=L+1;

end

end

附录二主程序

close all;clear all

s_f=15;s_xy=15;s_t=20;

s_mar=10;w_f=2.5;w_axi=2;

%信号产生

T0=8e-2;Ts=1/800;

N=T0/Ts;

t=0:Ts:(N-1)*Ts;

f1=50;f2=100;f3=200;fp=80;fs=120;

Sig=4*sin(f1*2*pi*t)+4*sin(f2*2*pi*t)+4*sin(f3*2*pi*t); s_need=4*sin(f1*2*pi*t);

s_1=4*sin(f2*2*pi*t);

s_2=4*sin(f3*2*pi*t);

P_qian=fftlzr(Sig);

figure(1)

h1=stem(2*pi/(N*Ts).*(0:N-1),abs(P_qian));grid on

hx1=xlabel('模拟角频率/ rad/s');

ht1=title('滤波前幅度频谱');

ha1=gca;

axis tight

%滤波器设计

Fs=1/Ts;

wp=2/Ts*tan(fp*2*pi*Ts/2);

ws=2/Ts*tan(fs*2*pi*Ts/2);

rp=1;rs=15;

[Nl,Wn]=buttord(wp,ws,rp,rs,'s');%求阶数及3dB截止频率[z,p,k]=buttap(Nl);%求极零点增益

[b,a]=zp2tf(z,p,k);%系统传递函数

[b2,a2]=lp2lp(b,a,Wn);%转换为所需低通滤波器

[HS,WS]=freqs(b2,a2,linspace(0,2000,512));%求0到2000rad/s的频率响应[hs,wss]=freqs(b2,a2,[wp,ws,[f1,f2,f3].*2*pi]);%求通带阻带边缘角频率及原始信号频率响应

%画模拟滤波器频谱并标注

figure(6)

subplot(211)

hsf=plot(WS,20*log10(abs(HS)));

hold on

stem(wss,20*log10(abs(hs)),'filled','r:','linewidth',w_f)

text(wss(1)-100,20*log10(abs(hs(1)))-10,...

{'\Omega_p=160\pi rad/s';[num2str(20*log10(abs(hs(1)))),'dB']}); text(wss(2)-65,20*log10(abs(hs(2)))-12,...

{'\Omega_s=240\pi rad/s';[num2str(20*log10(abs(hs(2)))),'dB']}); text(wss(3)-100,20*log10(abs(hs(3)))-10,...

{'\Omega_1=100\pi rad/s';[num2str(20*log10(abs(hs(3)))),'dB']}); text(wss(4)-65,20*log10(abs(hs(4)))-12,...

{'\Omega_2=200\pi rad/s';[num2str(20*log10(abs(hs(4)))),'dB']}); text(wss(5)-70,20*log10(abs(hs(5)))-10,...

{'\Omega_3=400\pi rad/s';[num2str(20*log10(abs(hs(5)))),'dB']}); hxsf=xlabel('模拟角频率/ rad/s');hysf=ylabel('幅值/dB');

htsf=title('幅频特性');hasf=gca;grid on

subplot(212)

hsx=plot(WS,angle(HS));

hold on

stem(wss,angle(hs),'filled','r:','linewidth',w_f)

text(wss(1),angle(hs(1))+0.7,...

{'\Omega_p=160\pi rad/s';[num2str(angle(hs(1))),'rad']});

text(wss(2),angle(hs(2))+0.7,...

{'\Omega_s=240\pi rad/s';[num2str(angle(hs(2))),'rad']});

text(wss(3)-150,angle(hs(3))-0.8,...

{'\Omega_1=100\pi rad/s';[num2str(angle(hs(3))),'rad']});

text(wss(4)-20,angle(hs(4))+0.8,...

{'\Omega_2=200\pi rad/s';[num2str(angle(hs(4))),'rad']});

text(wss(5)-50,angle(hs(5))-1.3,...

{'\Omega_3=400\pi rad/s';[num2str(angle(hs(5))),'rad']});

hxsx=xlabel('模拟角频率/ rad/s');hysx=ylabel('相值/ rad');

htsx=title('相频特性');hasx=gca;grid on

[b3,a3]=bilinear(b2,a2,Fs);%双线性变换到数字滤波器

[H,W]=freqz(b3,a3,linspace(0,0.5*pi,512));%求数字滤波器频率响应[hd,wd]=freqz(b3,a3,[wp,ws,[f1,f2,f3].*2*pi]*Ts);%求通带阻带边缘数字角频率及原始信号数字频率响应

wd=wd/pi;

%画数字滤波器频谱并标注

figure(2)

subplot(211)

h2=plot(W/pi,20*log10(abs(H)));

hold on

stem(wd,20*log10(abs(hd)),'filled','r:','linewidth',w_f) text(wd(1)-0.0035,20*log10(abs(hd(1)))-12,...

{['\omega_p=',num2str(wd(1)),'\pi

rad'];[num2str(20*log10(abs(hd(1)))),'dB']});

text(wd(2)-0.003,20*log10(abs(hd(2)))-12,...

{['\omega_s=',num2str(wd(2)),'\pi

rad'];[num2str(20*log10(abs(hd(2)))),'dB']});

text(wd(3)-0.0035,20*log10(abs(hd(3)))-12,...

{['\omega_1=',num2str(wd(3)),'\pi

rad'];[num2str(20*log10(abs(hd(3)))),'dB']});

text(wd(4)-0.003,20*log10(abs(hd(4)))-12,...

{['\omega_2=',num2str(wd(4)),'\pi

rad'];[num2str(20*log10(abs(hd(4)))),'dB']});

text(wd(5)-0.003,20*log10(abs(hd(5)))-12,...

{['\omega_3=',num2str(wd(5)),'\pi

rad'];[num2str(20*log10(abs(hd(5)))),'dB']});

hx2=xlabel('归一化数字角频率/ \pi rad');

hy1=ylabel('幅值/dB');

ht2=title('幅频特性');

ha2=gca;

grid on

subplot(212)

h3=plot(W/pi,angle(H));grid;title('相频特性'); hold on

stem(wd,angle(hd),'filled','r:','linewidth',w_f) text(wd(1),angle(hd(1))+1,...

{['\omega_p=',num2str(wd(1)),'\pi

rad'];[num2str(angle(hd(1))),'\pi rad']});

text(wd(2)-0.002,angle(hd(2))+1,...

{['\omega_p=',num2str(wd(2)),'\pi

rad'];[num2str(angle(hd(2))),'\pi rad']});

text(wd(3)-0.003,angle(hd(3))-1,...

{['\omega_p=',num2str(wd(3)),'\pi

rad'];[num2str(angle(hd(3))),'\pi rad']});

text(wd(4)-0.003,angle(hd(4))-2,...

{['\omega_p=',num2str(wd(4)),'\pi

rad'];[num2str(angle(hd(4))),'\pi rad']});

text(wd(5)-0.003,angle(hd(5))-1.3,...

{['\omega_p=',num2str(wd(5)),'\pi