matlab 快速傅立叶变换 FFT 实验报告Word格式.docx
- 文档编号:21164127
- 上传时间:2023-01-28
- 格式:DOCX
- 页数:8
- 大小:24.21KB
matlab 快速傅立叶变换 FFT 实验报告Word格式.docx
《matlab 快速傅立叶变换 FFT 实验报告Word格式.docx》由会员分享,可在线阅读,更多相关《matlab 快速傅立叶变换 FFT 实验报告Word格式.docx(8页珍藏版)》请在冰豆网上搜索。
k=1:
6,k=123456%fix向0舍入为整数,舍尾法x=cos(5*pi*n1*T1)+4*sin(2*pi*n1*T1).*sin(3*pi*n1*T1);
xa=cos(pi*n1*T1);
%obtainx(nT)andxa(nT)withgivensamplingfreq.'
fs'
fs=3;
T=1/fs;
n=0:
fix(time_period/T);
x_sample=cos(5*pi*n*T)+4*sin(2*pi*n*T).*sin(3*pi*n*T);
xa_sample=cos(pi*n*T);
figure,plot(n1*T1,x,'
r'
n1*T1,xa,'
b'
n*T,x_sample,'
ro'
),
%figure指令可以打开一个新的图形绘制窗口,在下一条这样的指令出现之前,
%所有的绘图指令都会在这个窗口上绘制。
%plot函数可以同时画多组数值图像,即以(X坐标数组1,Y坐标数组1,X%坐标数组2,Y坐标数组2,……)的形式书写。
每组坐标书写后,可用一个‘’
%表示前面的一组坐标用怎样的形式绘制:
线型,点型和颜色。
如‘r’就表示%用绿色的线绘制,‘b’蓝色。
‘ro’红色圆圈等等。
holdon,stem(n*T,xa_sample,'
b:
x'
)%hold保持当前图形窗的内容即:
在画完一张图后用hold命令保持住,继续在%当前窗画图。
%stem离散序列绘图%b:
x表示蓝色虚线x号legend('
x(t)'
'
xa(t)'
x(nT)'
xa(nT)'
),xlabel('
t(ms)'
)%legend这个指令可以在图像的右上角绘出一个图例,表示每条曲线代表什么,%名称需编程者指定。
%xlabelx轴标注,ylabel:
y轴标注,title:
三维坐标标注。
上面程序结果如下图所示:
波形分析:
x(t)与xa(t)两波形并不相同,但采样后,波形一致。
二、判别离散时间系统的时不变性。
p105例3.2.2)
设输入序列为EMBEDEquation.3,系统EMBEDEquation.3实现对EMBEDEquation.3的抽取。
设EMBEDEquation.3。
取延迟量D(例如D=30)。
记EMBEDEquation.3,画出EMBEDEquation.3、EMBEDEquation.3的序列波形。
编程求出系统对EMBEDEquation.3的响应EMBEDEquation.3以及对EMBEDEquation.3的响应EMBEDEquation.3画出EMBEDEquation.3、EMBEDEquation.3的波形。
该系统是否为时不变的?
%=============%problem2%=============clear%plotx(n)andx(n-D)D=30;
N=500;
n=1:
N;
%增量为1的等增量语句x=sin(2*pi/100*n);
forn=1:
N+D,if(n-D)<
=0,xD(n)=0;
elsexD(n)=x(n-D);
endend%求xD(n)的表达式figure,subplot(2,1,1),%subplot(n,m,p),将图形分成n*m个子图,在第p个子图内绘图plot(1:
N,x,'
r:
'
1:
length(xD),xD,'
),legend('
x(n)'
xD(n)'
n'
)%ploty(n)andyD(n)forn=1:
fix(N/2)%fix,对N/2去尾取整y(n)=x(2*n);
endforn=1:
length(y)+D,if(n-D)<
=0,y_delay(n)=0;
elsey_delay(n)=y(n-D);
endend%输入x(n),输出延迟Dforn=1:
fix(length(xD)/2)yD(n)=xD(2*n);
end%输入延迟D,对应的此时的输出subplot(2,1,2),%在第二个子图内绘图plot(1:
length(y),y,'
length(y_delay),y_delay,'
r.:
length(yD),yD,'
b.'
y(n)'
y(n-D)'
yD(n)'
),,xlabel('
)axis([0530-11])%axis(v),v四元向量[xmin,xmax,ymin,ymax],坐标轴设定在v规定的范围内上面程序结果如下图所示:
根据系统时不变的性质,当xD(n)=x(n-D)时,有yD(n)=y(n-D)。
由图2可以看出,y(n-D)与yD(n)的波形并不一致,则:
系统为时变系统
三、利用卷积计算信号通过FIR滤波器的输出,并观察输出信号的input-on暂态、input-off暂态和稳态阶段。
p144例4.1.8)
考虑两个滤波器,EMBEDEquation.3,EMBEDEquation.3;
输入EMBEDEquation.3为周期方波,第一个周期内EMBEDEquation.3。
1..分别画出EMBEDEquation.3通过两个滤波器的输出EMBEDEquation.3、EMBEDEquation.3的波形EMBEDEquation.3,并与书上p144例4.1.8的两幅图比较是否一致。
2.计算图中稳态部分的响应值。
%=============%problem3%=============clearh1=0.25*0.75.^(0:
14);
h2=1/5*[1-510-105-1];
%矩阵表示方法N=200;
n=0:
N-1;
%增量为1的等增量语句x1=[ones(1,25)zeros(1,25)];
%oneperiodof'
%ones(m,n),全1矩阵(m行,n列);
zeros(m,n),全0矩阵(m行,n列)%eye(n),单位矩阵,n阶方阵x=[x1x1x1x1];
y1=conv(x,h1);
y2=conv(x,h2);
%y=conv(u,h)卷积语句figure,subplot(2,1,1),,plot(n,x,'
n,y1(1:
N),'
),%将图形分成2个子图,在第一个子图内绘图axis([0200-0.52.5]),gridon,%grid:
图上加坐标网格legend('
input'
output'
)subplot(2,1,2),plot(n,x,'
n,y2(1:
),axis([0200-1.52.5]),gridon,%图上加坐标网格legend('
)上面程序结果如下图所示:
与书上P144图形一致2.1EMBEDEquation.3
ItsDCgainisalmostunity:
(M=14,m<
=M且m>
=0)
Ydc=∑h(m)=0.25∑0.75^m=1-0.75^(14+1)=0.987
2.2EMBEDEquation.3H(Z)=1/5(1-1/Z)^5ydc=H(Z)|z=1→ydc=0
实验二
一、观察序列频谱,观察信号通过系统后波形与频谱的变化
已知输入信号EMBEDEquation.3,其中EMBEDEquation.3,EMBEDEquation.3,EMBEDEquation.3,N可取5000点。
(1)画出EMBEDEquation.3的前100点波形
(2)画出EMBEDEquation.3的DTFT频谱EMBEDEquation.3(EMBEDEquation.3)
由于计算机无法画出连续频谱,所以可在EMBEDEquation.3内均匀取足够密的点数,如M=5000个频率点EMBEDEquation.3EMBEDEquation.3,求出这些频率点上的频谱值EMBEDEquation.3EMBEDEquation.3,并画出EMBEDEquation.3随EMBEDEquation.3变化的曲线。
(3)某LTI系统EMBEDEquation.3,画出系统的幅度频响EMBEDEquation.3
(4)求系统对EMBEDEquation.3的响应EMBEDEquation.3(可以自己编程也可利用卷积函数EMBEDEquation.3)。
画出EMBEDEquation.3的波形,并与EMBEDEquation.3的波形比较(各画100点);
画出EMBEDEquation.3的幅度谱EMBEDEquation.3,并与EMBEDEquation.3比较。
你从中观察到什么?
Promble1
二、系统函数EMBEDEquation.3,根据正准型结构(canonicalform)编写样本处理算法。
内部状态的初始值设为零,输入信号EMBEDEquation.3采用逐个样本手动输入的方式,求输出信号EMBEDEquation.3。
Promble2w1=0;
fori=1:
10,x=input('
inputx='
);
w0=0.8*w1+x;
y=w0-0.5*w1;
w1=w0;
yendinputx=0y=0inputx=1y=1inputx=2y=2.3000inputx=3y=3.8400inputx=4y=5.5720inputx=5y=7.4576inputx=6y=9.4661inputx=7y=11.5729inputx=8y=13.7583inputx=9y=16.0066
实验三
理解例8.3.1(p393)IIR平滑器的作用和特点(来源:
p468习题8.30)
(1)重现p394图8.3.4中a=0.90和a=0.98两情形下L=200点的输出响应;
EMBEDEquation.3,其中s=5,噪声方差EMBEDEquation.3
相关程序:
%--------Problem1---------------clearL=200;
%(用来确定信号的长度)s=5*ones(L,1);
%(ones()这个函数用来产生单位序列,高度为5)v=randn(L,1);
%(randn()是用来产生随机序列,喜爱MATLAB中这个函数是以均值为零方差为一的方式产生随机序列。
)v=v-mean(v);
%(下面这两句话是对产生的长度为200的序列进行修正,因为随机序列v虽然是以高斯白噪声产生的,不过由于其随机性,有可能会造成均值和方差不是0,1所以要分别减去均值mean(v),方差std(v))v=v/std(v);
x=s+v;
%(此时有用信号加上模拟出来的白噪声信号v,即是信道传输出来的信号)%(下面一步就是对滤波器的模拟,注意filter的用法:
filter(b,【1-a】,x)这是一种差分方程的表示方法y(n)-a*y(n-1)=b*x(n)用这种方法可以比较方便的表示一个差分方程及其响应)a=0.90;
b=1-a;
y1=filter(b,[1-a],x)'
;
%(注意!
方框中不是1-a而是1-a)%(接下来就是对a=0.98的响应y2,并且进行画图)a=0.98;
y2=filter(b,[1-a],x)'
subplot(2,1,1),plot((1:
L),x,'
:
(1:
L),y1,(1:
L),s,'
),title('
a=0.90'
)subplot(2,1,2),plot((1:
L),y2,(1:
legend('
a=0.98'
)
Ps:
图中的噪声信号和老师的信号有所不同原因是randn产生的噪声信号不同,不过总体趋势是一样的!
(2)估算NRR。
估算公式为:
输入噪声方差估值:
EMBEDEquation.3,输出噪声方差估值:
EMBEDEquation.3
NRR估值:
mx=mean(x);
varx=sum((x-mx).^2)/L;
%(请大家注意:
(x-mx)后面有‘.’是对矩阵运算时才愈要加上!
根据公式估算出方差)start=1;
%(从y=1这一点开始计算方差,会发现计算值与理论值有很大的误差)a=0.90;
y1=y1(start:
L);
my=mean(y1);
vary=sum((y1-my).^2)/length(y1);
NRResti=vary/varx;
NRR=(1-a)/(1+a);
[aNRRestiNRR]a=0.98;
y2=y2(start:
my=mean(y2);
vary=sum((y2-my).^2)/length(y1);
[aNRRestiNRR]
根据以上的第一问的噪声信号所计算的理论值和估算值:
(3)计算NRR理论值,与估计值比较。
为何相差较大?
如果想估计得准确,该怎么办?
(提示:
增大L,增大程序中变量start)
Start=100时:
Start=100;
L=2000时:
可以看出增大start可以减少与L的差别,而增大L则会使NRR减小!
2.(来源p469习题8.31)
仍用上述一阶IIR平滑器,只是输入高频信号EMBEDEquation.3
b取何值时,信号部分通过系统后不变(指稳态)?
如果你已经完成了作业,验证一下。
代码如下:
%----------Problem2--------------clearL=200;
s=5*(-1).^(1:
%cleansignalv=randn(1,L);
v=v-mean(v);
v=v/std(v);
%(与上题一样v是一个高斯白噪声信号)x=s;
%(首先输入一个无杂波的信号,看系统响应,是否能够产生出输入信号)a=0.98;
b=1+a;
y=filter(b,[1-a],x);
figure,plot((1:
L),y),
s(n)'
),axis([1L-106]),title('
nonoise'
从此可以看出,系统的冲击响应在与输入信号S(n)卷积后经过一段时间能还原出输入信号,即此滤波器可以实现原始信号的通过。
(2)噪声部分是否放大了?
然后再把噪声加进去:
x=s+v;
a=0.98;
L),y),legend('
)
这时,我们看到通过滤波器的信号已经无法恢复出输入信号了,证明滤波器无法抑制噪声信号,
我们在进行编程看系统是否对噪声有放大作用:
NRR=b.^2/(1-a.^2)
NRR=99.0000
系统对信号放大了99倍!
!
M=2500;
w=pi/M*(0:
M-1);
j=sqrt(-1);
k=1:
M,Hw=b./(1-a*exp(-j*w(k)));
%(求出H(w)即使用公式EMBEDEquation.KSEE3\*MERGEFORMAT从而得到系统的频谱)H=abs(Hw).^2;
figure,plot(w/pi,H,w/pi,ones(M,1),'
|H(w)|^2'
inputnoisespectrum'
\omega(\pi)'
),ylabel('
|H(\omega)|^2'
小结:
通过这次试验,和自己的对照测试,查‘help’觉得收获最大的就是了解了filter()这个函数的用法。
很方便的根据差分方程或者z域里的方程直接编程写出,而且效率很高!
实验四
观察窗函数的影响。
信号为EMBEDEquation.3,EMBEDEquation.3KHz,EMBEDEquation.3KHz,EMBEDEquation.3KHz,采样频率EMBEDEquation.3KHz。
写出EMBEDEquation.3(EMBEDEquation.3)的频谱EMBEDEquation.3;
分别画出窗长度EMBEDEquation.3,EMBEDEquation.3,EMBEDEquation.3,EMBEDEquation.3的矩形窗频谱和Hamming窗频谱。
观察L的变化对窗谱的主瓣宽度、旁瓣密集度、相对旁瓣水平的影响。
c)时域采样点数分别取EMBEDEquation.3,EMBEDEquation.3,EMBEDEquation.3,EMBEDEquation.3,
分别画出EMBEDEquation.3加矩形窗及加Hamming窗时DTFT频谱EMBEDEquation.3;
你能否从频谱上分辨出信号的三个频率分量?
若能分辨出,它们的位置和相对大小是否准确?
L的大小对频率的物理分辨率(physicalfrequencyresolution)有何影响?
a):
clear;
N=4000;
fs=10000;
t=0:
1/fs:
(N-1)/fs;
w=0:
2*pi/N:
(N-1)/N*2*pi;
f1=2000;
f2=2500;
f3=3000;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t)+cos(2*pi*f3*t);
X=fft(x);
plot(w/pi,abs(X)/N);
xlabel('
\omegainunitsof\pi'
ylabel('
|X(\omega)|'
L=10;
subplot(2,2,1);
windows_spectrum(L);
L=20;
subplot(2,2,2);
L=40;
subplot(2,2,3);
L=100;
subplot(2,2,4);
L=100;
x=cos(0.4*pi*(0:
(L-1)))+cos(0.5*pi*(0:
(L-1)))+cos(0.6*pi*(0:
(L-1)));
windowed_spectrum(x(1:
L));
由图形可看出,当L40时,汉明窗才可以分辨出三个峰值,而矩形窗只要L20即可。
理解频率的物理分辨率和计算分辨率的区别
信号同上,加矩形窗。
时域采样点数分别取EMBEDEquation.3,EMBEDEquation.3,EMBEDEquation.3,EMBEDEquation.3。
画出以上各种时长情况下,频域采样点数分别为EMBEDEquation.3,EMBEDEquation.3时的DFT(在同一个图上用虚线画出相应的DTFT频谱,用于比较)。
离散频谱DFT和连续频谱DTFT有什么关系?
L一定的情况下,能否通过增加N改善频率的物理分辨率?
N的作用是什么?
N0=32;
windowed_dft_spectrum(x(1:
L),N0);
DFT相当于DTFT在频域的采样。
将两张图对比可发现,当L较小(如10),此时N增加也无效果,因为物理分辨率由L决定。
N决定在频域采样的密度,也就是决定了计算分辨率。
N0=64;
另附三个函数的源代码。
由于时间所限,我减少了老师提供的程序之中的输入参数,使得这几个函数只提供对汉明窗和矩形窗的分析。
以及附上的功能函数的有关解释:
functionwindows_spect
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 快速傅立叶变换 FFT 实验报告 快速 傅立叶 变换 实验 报告