matlab报告Word文件下载.docx
《matlab报告Word文件下载.docx》由会员分享，可在线阅读，更多相关《matlab报告Word文件下载.docx（20页珍藏版）》请在冰豆网上搜索。
通过c语言实现的快速傅里叶变换具有较为复杂的编程过程，这对与计算具有很大的不便，但是如果采用matlab的话就比较简单了。
第四次实验：
IIR和FIR滤波器的matlab实现
通过matlab实现特定要求的滤波器，例如巴特沃斯、切比雪夫，还有FIR滤波器的实现
切比雪夫带通滤波器
巴特沃斯带通滤波器
通过比较就可以得出互相的优缺点，而且过程较为简单。
1.滤波器的幅频与相频响应
2.输入与滤波后信号显示
3.输出的信号波形
各种滤波窗口的比较：
实验总结：
通过计算，编程可以得出滤波器的滤波特性。
附录：
三个同心圆的显示：
t=0:
0.01:
2*pi;
x=exp（i*t）;
y=[x;
2*x;
3*x]'
;
plot（y）
gridon;
%加网格线
boxon;
%加坐标边框
axisequal%坐标轴采用等刻度
饼状图比例分析：
subplot（1,2,1）;
pie（[2347,1827,2043,3025]）;
title（'
饼图'
）;
legend（'
一季度'
'
二季度'
三季度'
四季度'
subplot（1,2,2）;
compass（[7+2.9i,2-3i,-1.5-6i]）;
相量图'
动态图形显示：
[X,Y,Z]=peaks（30）;
surf（X,Y,Z）
axis（[-3,3,-3,3,-10,10]）
axisoff;
shadinginterp;
colormap（hot）;
m=moviein（20）;
%建立一个20列大矩阵
fori=1:
20
view（-37.5+24*（i-1）,30）%改变视点
m（:
i）=getframe;
%将图形保存到m矩阵
end
movie（m,2）;
%播放画面2次
卷积的计算公式：
function[y,ny]=conv_m（x,nx,h,nh）
nyb=nx
（1）+nh
（1）;
nye=nx（length（x））+nh（length（h））;
ny=[nyb:
nye];
y=conv（x,h）;
卷积的应用实例：
x=[00.511.50];
nx=0:
4;
h=[11100];
nh=0:
[y,ny]=conv_m（x,nx,h,nh）;
subplot（2,2,1）;
stem（nx,x）;
序列x'
xlabel（'
n'
ylabel（'
x（n）'
subplot（2,2,2）;
stem（nh,h）;
序列h'
h（n）'
subplot（2,2,3）;
stem（ny,y）;
两序列卷积'
y（n）'
巴特沃斯滤波器的实现
w1=2*400*tan（2*pi*95/（2*400））;
w2=2*400*tan（2*pi*105/（2*400））;
wr=2*400*tan（2*pi*110/（2*400））;
wm=2*400*tan（2*pi*90/（2*400））;
[n,wn]=buttord（[w1w2],[wmwr],3,10,'
s'
[b,a]=butter（n,wn,'
[num,den]=bilinear（b,a,400）;
[h,w]=freqz（num,den）;
f=w/pi*200;
plot（f,20*log10（abs（h）））;
axis（[90,110,-10,3]）;
grid;
频率/kHz'
）
幅度/dB'
切比雪夫滤波器的实现
[n,wn]=cheb1ord（[w1w2],[wmwr],3,10,'
[b,a]=cheby1（n,3,wn,'
双线性变换实现巴特沃斯低通滤波器的技术指标：
采样频率10Hz。
通带截止频率fp=0.2*piHz。
阻带截止频率fs=0.3*piHz。
4.通带衰减小于1dB，阻带衰减大于20dB
使用双线性变换法由模拟滤波器原型设计数字滤波器
T=0.1;
FS=1/T;
fp=0.2*pi;
fs=0.3*pi;
wp=fp/FS*2*pi;
ws=fs/FS*2*pi;
Rp=1;
%通带衰减
As=15;
%阻带衰减
OmegaP=（2/T）*tan（wp/2）;
%频率预计
OmegaS=（2/T）*tan（ws/2）;
%设计巴特沃斯低通滤波器原型
N=ceil（（log10（（10^（Rp/10）-1）/（10^（As/10）-1）））/（2*log10（OmegaP/OmegaS）））;
OmegaC=OmegaP/（（10^（Rp/10）-1）^（1/（2*N）））;
[z,p,k]=buttap（N）;
%获取零极点参数
p=p*OmegaC;
k=k*OmegaC^N;
B=real（poly（z））;
b0=k;
cs=k*B;
ds=real（poly（p））;
[b,a]=bilinear（cs,ds,FS）;
%双线性变换
figure
（1）;
%绘制结果
freqz（b,a,512,FS）;
%进行滤波验证
figure
（2）;
%绘制结果
figure（3）;
f1=50;
f2=250;
n=0:
63;
x=sin（2*pi*f1*n）+sin（2*pi*f2*n）;
stem（x,'
.'
title（'
输入信号'
y=filter（b,a,x）;
stem（y,'
）;
title（'
滤波之后的信号'
输出的信号'
wc=0.2*pi;
N=33;
tao=（N-1）/2;
n=[0:
（N-1）];
m=n-tao+0.000001;
hd=sin（wc*m）./（pi*m）;
wd1=boxcar（N）'
b1=hd.*wd1;
wd2=hanning（N）'
b2=hd.*wd2;
wd3=blackman（N）'
b3=hd.*wd3;
wd4=hamming（N）'
b4=hd.*wd4;
[h1,w]=freqz（b1,1）;
[h2,w]=freqz（b2,1）;
[h3,w]=freqz（b3,1）;
[h4,w]=freqz（b4,1）;
plot（w,20*log10（abs（h1））,’-’,20*log10（abs（h2））,’..’,20*log10（abs（h3））,’-.’,20*log10（abs（h4））,’--‘）;
矩形窗'
汉宁窗'
布莱克曼窗'
汉明窗'
c语言实现fft
#include<
stdio.h>
dos.h>
stdlib.h>
malloc.h>
math.h>
typedefstruct{
floatreal;
floatimag;
}COMPLEX;
externvoidfft（COMPLEX*x,intm）;
voidfft（COMPLEX*x,intm）
{
staticCOMPLEX*w;
staticintmstore=0;
staticintn=1;
COMPLEXu,temp,tm;
COMPLEX*xi,*xip,*xj,*wptr;
inti,j,k,l,le,windex;
doublearg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real,wtemp_imag;
if（m!
=mstore）
{
if（mstore!
=0）free（w）;
mstore=m;
if（m==0）return;
n=1<
<
m;
le=n/2;
w=（COMPLEX*）malloc（（le-1）*sizeof（COMPLEX））;
if（!
w）
{
printf（"
\nUnabletoallocateCOMPLEXwarray\n"
exit
（1）;
}
arg=4.0*atan（1.0）/le;
wrecur_real=w_real=cos（arg）;
wrecur_imag=w_imag=-sin（arg）;
xj=w;
for（j=1;
j<
le;
j++）
xj->
real=（float）wrecur_real;
imag=（float）wrecur_imag;
xj++;
wtemp_real=wrecur_real*w_real-wrecur_imag*w_imag;
wtemp_imag=wrecur_real*w_imag+wrecur_imag*w_real;
wrecur_real=wtemp_real;
}
le=n;
windex=1;
for（l=0;
l<
l++）
le=le/2;
for（i=0;
i<
n;
i=i+2*le）
xi=x+i;
xip=xi+le;
temp.real=（xi->
real+xip->
real）;
temp.imag=（xi->
imag+xip->
imag）;
xip->
real=（xi->
real-xip->
imag=（xi->
imag-xip->
xi=&
temp;
wptr=w+windex-1;
u=*wptr;
for（i=j;
{
xi=x+i;
xj=x+j;
temp=*xj;
*xj=*xi;
*xi=temp;
}
}
voidmain（）
COMPLEX*x,*xtemp;
inti,nn,mm;
doubley,q;
system（"
cls"
printf（"
\nInputm="
scanf（"
%d"
&
mm）;
nn=1<
mm;
x=（COMPLEX*）malloc（nn*sizeof（COMPLEX））;
xtemp=（COMPLEX*）malloc（nn*sizeof（COMPLEX））;
x[0].real=1.0;
x[1].imag=0.0;
for（i=1;
nn;
i++）
x[i].real=（float）exp（-i）;
x[i].imag=0.0;
for（i=0;
i++）xtemp[i]=x[i];
fft（x,mm）;
Finished!
\n\n"
kXR（K）XI（K）|X（K）|Q（k）\n"
y=x[i].real*x[i].real+x[i].imag*x[i].imag;
if（x[i].real!
=0）
if（x[i].imag==0.0）
if（x[i].real>
0）q=0.0;
if（x[i].real<
0）q=4.0*（atan（1.0））;
q=atan（x[i].imag/x[i].real）;
printf（"
%4d%9f%9f%9f%9f\n"
i,x[i].real,x[i].imag,sqrt（y）,q）;
通过四次实验，我们并没有完成太多的任务，只是简单的验证了一下书上的一些代码，我认为这门课程还是很有必要多学习一下的，这对于以后进行有关系统建立，信号处理，等都有十分大的帮助，所以我们在平时也应该多学习有关该软件的有关使用。
升级会员 