收藏
下载资源下载资源 加入VIP,免费下载

信号分析与处理实验报告.docx

上传人：b****6 文档编号：8252879 上传时间：2023-01-30 格式：DOCX 页数：20 大小：217.55KB
下载 相关 举报
信号分析与处理实验报告.docx_第1页
第1页 / 共20页
信号分析与处理实验报告.docx_第2页
第2页 / 共20页
信号分析与处理实验报告.docx_第3页
第3页 / 共20页
信号分析与处理实验报告.docx_第4页
第4页 / 共20页
信号分析与处理实验报告.docx_第5页
第5页 / 共20页
点击查看更多>>
下载资源
资源描述

信号分析与处理实验报告.docx

《信号分析与处理实验报告.docx》由会员分享，可在线阅读，更多相关《信号分析与处理实验报告.docx（20页珍藏版）》请在冰豆网上搜索。

信号分析与处理实验报告.docx

信号分析与处理实验报告

 

华北电力大学

实验报告

|

|

 

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

IIR数字滤波器的设计

课程名称信号分析与处理

|

|

专业班级：

电气化1308学生姓名：

袁拉麻加

学号：

********0827成绩：

指导教师：

杨光实验日期：

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

程序：

%%

clearall;closeall;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

%%

%码位倒置%

fork=0:

n-1

forj=1:

m%取M位的二进制数%

x1（j）=bitget（k,j）;%倒取出二进制数%

end

x1=num2str（x1）;%将数字序列转化为字符串%

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

clearx1

end

fork=1:

n

B（k）=A（y（k）+1）;%时间抽取序列%

end

clearA

%%

%计算%

forL=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）;

forR=1:

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

forP=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

程序：

%%

clearall;closeall;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

%%

%码位倒置%

fork=0:

n-1

forj=1:

m%取M位的二进制数%

x1（j）=bitget（k,j）;%倒取出二进制数%

end

x1=num2str（x1）;%将数字序列转化为字符串%

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

clearx1

end

fork=1:

n

B（k）=A（y（k）+1）;%时间抽取序列%

end

clearA

%%

%计算%

forL=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）;

forR=1:

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

forP=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的过程更加熟悉了，尤其是WN系数的特性，可以通过累乘来实现W系数的更新。

IIR数字滤波器的设计实验

一、实验目的

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

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

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

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

二、实验内容

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

、

、

，将它们叠加成一个信号

。

用满足采样定理的采样频率对模拟信号

进行采样，得到离散时间序列

。

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

进行FFT运算，观察它的频谱。

3.设计满足已知的设计指标的IIR数字滤波器，并用它对

进行滤波处理。

对滤波后得到的波形再做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）;

iffix（M）~=M

A=[A（:

）;zeros（2^（fix（M）+1）-N,1）];

end

N=length（A）;M=log2（N）;

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

whileI~=N-1

ifI

t=A（J+1）;

A（J+1）=A（I+1）;

A（I+1）=t;

end

K=N/2;

whileK<=J

J=J-K;

K=K/2;

end

J=J+K;

I=I+1;

end

whileL<=M

LE=2^L;R=1;

LE1=LE/2;

W=1;W1=exp（-j*pi/LE1）;

whileR<=LE1

P=R;

whileP<=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

附录二主程序

closeall;clearall

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））;gridon

hx1=xlabel（'模拟角频率/rad/s'）;

ht1=title（'滤波前幅度频谱'）;

ha1=gca;

axistight

%滤波器设计

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）））;

holdon

stem（wss,20*log10（abs（hs））,'filled','r:

','linewidth',w_f）

text（wss

（1）-100,20*log10（abs（hs

（1）））-10,...

{'\Omega_p=160\pirad/s';[num2str（20*log10（abs（hs

（1））））,'dB']}）;

text（wss

（2）-65,20*log10（abs（hs

（2）））-12,...

{'\Omega_s=240\pirad/s';[num2str（20*log10（abs（hs

（2））））,'dB']}）;

text（wss（3）-100,20*log10（abs（hs（3）））-10,...

{'\Omega_1=100\pirad/s';[num2str（20*log10（abs（hs（3））））,'dB']}）;

text（wss（4）-65,20*log10（abs（hs（4）））-12,...

{'\Omega_2=200\pirad/s';[num2str（20*log10（abs（hs（4））））,'dB']}）;

text（wss（5）-70,20*log10（abs（hs（5）））-10,...

{'\Omega_3=400\pirad/s';[num2str（20*log10（abs（hs（5））））,'dB']}）;

hxsf=xlabel（'模拟角频率/rad/s'）;hysf=ylabel（'幅值/dB'）;

htsf=title（'幅频特性'）;hasf=gca;gridon

subplot（212）

hsx=plot（WS,angle（HS））;

holdon

stem（wss,angle（hs）,'filled','r:

','linewidth',w_f）

text（wss

（1）,angle（hs

（1））+0.7,...

{'\Omega_p=160\pirad/s';[num2str（angle（hs

（1）））,'rad']}）;

text（wss

（2）,angle（hs

（2））+0.7,...

{'\Omega_s=240\pirad/s';[num2str（angle（hs

（2）））,'rad']}）;

text（wss（3）-150,angle（hs（3））-0.8,...

{'\Omega_1=100\pirad/s';[num2str（angle（hs（3）））,'rad']}）;

text（wss（4）-20,angle（hs（4））+0.8,...

{'\Omega_2=200\pirad/s';[num2str（angle（hs（4）））,'rad']}）;

text（wss（5）-50,angle（hs（5））-1.3,...

{'\Omega_3=400\pirad/s';[num2str（angle（hs（5）））,'rad']}）;

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

htsx=title（'相频特性'）;hasx=gca;gridon

[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）））;

holdon

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））,'\pirad'];[num2str（20*log10（abs（hd

（1））））,'dB']}）;

text（wd

（2）-0.003,20*log10（abs（hd

（2）））-12,...

{['\omega_s=',num2str（wd

（2））,'\pirad'];[num2str（20*log10（abs（hd

（2））））,'dB']}）;

text（wd（3）-0.0035,20*log10（abs（hd（3）））-12,...

{['\omega_1=',num2str（wd（3））,'\pirad'];[num2str（20*log10（abs（hd（3））））,'dB']}）;

text（wd（4）-0.003,20*log10（abs（hd（4）））-12,...

{['\omega_2=',num2str（wd（4））,'\pirad'];[num2str（20*log10（abs（hd（4））））,'dB']}）;

text（wd（5）-0.003,20*log10（abs（hd（5）））-12,...

{['\omega_3=',num2str（wd（5））,'\pirad'];[num2str（20*log10（abs（hd（5））））,'dB']}）;

hx2=xlabel（'归一化数字角频率/\pirad'）;

hy1=ylabel（'幅值/dB'）;

ht2=title（'幅频特性'）;

ha2=gca;

gridon

subplot（212）

h3=plot（W/pi,angle（H））;grid;title（'相频特性'）;

holdon

stem（wd,angle（hd）,'filled','r:

','linewidth',w_f）

text（wd

（1）,angle（hd

（1））+1,...

{['\omega_p=',num2str（wd

（1））,'\pirad'];[num2str（angle（hd

（1）））,'\pirad']}）;

text（wd

（2）-0.002,angle（hd

（2））+1,...

{['\omega_p=',num2str（wd

（2））,'\pirad'];[num2str（angle（hd

（2）））,'\pirad']}）;

text（wd（3）-0.003,angle（hd（3））-1,...

{['\omega_p=',num2str（wd（3））,'\pirad'];[num2str（angle（hd（3）））,'\pirad']}）;

text（wd（4）-0.003,angle（hd（4））-2,...

{['\omega_p=',num2str（wd（4））,'\pirad'];[num2str（angle（hd（4）））,'\pirad']}）;

text（wd（5）-0.003,angle（hd（5））-1.3,...

{['\omega_p=',num2str（wd（5））,'\pirad'];[num2str（angle（hd（5）））,'\pirad']}）;

hx3=xlabel（'归一化数字角频率/\pirad'）;

hy2=ylabel（'相值/\pirad'）;

ht3=title（'相频特性'）;

ha3=gca;

gridon

%滤波

data=filter（b3,a3,Sig）;

figure（3）

h4=plot（t,data,'k',t,Sig,'b',t,s_need,'r'）;

ht4=title（'滤波效果'）;

xlabel（'t/s','fontsize',s_xy）;

ylabel（'幅值','fontsize',s_xy）;

legend（'滤波后波形','滤波前波形','50Hz信号'）

ha4=gca;

gridon

P_hou=fftlzr（data）;

figure（4）

subplot（211）

h5=stem（2*pi/（N*Ts）.*（0:

N-1）,abs（P_qian））;

hx4=xlabel（'模拟角频率/rad/s'）;

ht5=title（'滤波前幅度频谱'）;

axis（[0,1.5e3,0,2e3]）

text（f1*2*pi-20,abs（P_qian（f1*N*Ts+1））+200,[num2str（f1）'Hz']）;

text（f2*2*pi-20,abs（P_qian（f2*N*Ts+1））+200,[num2str（f2）'Hz']）;

text（f3*2*pi-20,abs（P_qian（f3*N*Ts+1））+200,[num2str（f3）'Hz']）;

ha5=gca;gridon

subplot（212）

h6=stem（2*pi/（N*Ts）.*（0:

N-1）,abs（P_hou））;

hx5=xlabel（'模拟角频率/rad/s'）;

ht6=title（'滤波后幅度频谱'）;

axis（[0,1.5e3,0,2e3]）

text（f1*2*pi-20,abs（P_hou（f1*N*Ts+1））+200,[num2str（f1）'Hz']）;

text（f2*2*pi-20,abs（P_hou（f2*N*Ts+1））+200,[num2str（f2）'Hz']）;

text（f3*2*pi-20,abs（P_hou（f3*N*Ts+1））+200,[num2str（f3）'Hz']）;

ha6=gca;gridon

figure（5）

h7=plot（t,Sig,t,s_need,'r',t,s_1,'g',t,s_2,'c'）;

hx6=xlabel（'t/s'）;

hy3=xlabel（'幅值'）;

ht7=title（'生成信号'）;

set（h7,'linewidth',2）

legend（'合成信号','50Hz','100Hz','200Hz'）;

ha7=gca;gridon

%图像处理

h=[h1,h2,h3,h4',h5,h6,hsf,hsx];

ht=[ht1,ht2,ht3,ht4,ht5,ht6,ht7,htsf,htsx];

hx=[hx1,hx2,hx3,hx4,hx5,hx6,hxsf,hxsx];

hy=[hy1,hy2,hy3,hysf,hysx];ha=[ha1,ha2,ha3,ha4,ha5,ha6,ha7,hasf,hasx];

set（h（[2:

6,9:

10]）,'Linewidth',w_f）

set（ha,'Linewidth',w_axi,'Fontsize',s_f）;

set（[hxhyht],'Fontsize',s_xy,'FontWeight','bold'）;

set（（1:

6）,'outerposition',get（0,'screensize'）,'color','w'）

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置：首页 > 小学教育 > 语文

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1