matlab实现DFT.docx
《matlab实现DFT.docx》由会员分享，可在线阅读，更多相关《matlab实现DFT.docx（14页珍藏版）》请在冰豆网上搜索。
matlab实现DFT
DFT基于Matlab的实现
一、实验目的
1．掌握DFT函数的用法。
2.利用DFT进行信号检测及谱分析。
3．了解信号截取长度对谱分析的影响。
二、实验内容
1．利用DFT计算信号功率谱。
实验程序：
t=0:
0.001:
0.6;
x=sin（2*pi*50*t）+sin（2*pi*120*t）+randn（1,length（t））;
Y=dft（x,512）;
P=Y.*conj（Y）/512;
f=1000*（0:
255）/512;
plot（f,P（1:
256））
2.进行信号检测。
分析信号频谱所对应频率轴的数字频率和频率之间的关系。
模拟信号
，以
进行取样，求N点DFT的幅值谱。
实验程序：
subplot（2,2,1）
N=45;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）;
y=dft（x,N）;
plot（q,abs（y））;title（'DFTN=45'）
subplot（2,2,2）
N=50;n=0:
N-1;t=0.01*n;q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）;
y=dft（x,N）;
plot（q,abs（y））;title（'DFTN=50'）
subplot（2,2,3）
N=55;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）;
y=dft（x,N）;
plot（q,abs（y））;title（'DFTN=55'）
subplot（2,2,4）
N=60;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）;
y=dft（x,N）;
plot（q,abs（y））;title（'DFTN=60'）
3.对2，进一步增加截取长度和DFT点数，如N加大到256，观察信号频谱的变化，分析产生这一变化的原因。
在截取长度不变的条件下改变采样频率，观察信号频谱的变化，分析产生这一变化的原因。
N加大到256时的程序：
N=256;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）;
y=dft（x,N）;
plot（q,abs（y））;title（'DFTN=256'）
分析原因：
在T=0.01s的情况下，第一个序列的周期是100，第二个序列的周期是50，所以当取样点数小于100时，频率分辨率不够，不能够区分出两个信号。
当采样点数足够多（256）时，频率分辨率增加，能够区分出两个频率的信号。
将采样间隔变为T=0.1s时，N仍为45的程序：
N=45;n=0:
N-1;t=0.1*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）;
y=dft（x,N）;
plot（q,abs（y））;title（'DFTN=45'）
分析原因：
在T=0.1s的情况下，第一个序列的周期是10，第二个序列的周期是5，所以当取样点数为45时，能够区分出两个信号。
参数同上，N取64，并在信号中加入噪声w（t）。
figure
（2）
subplot（2,1,1）
N=64;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）;y=dft（x,N）;
plot（q,abs（y））;title（'DFTN=64'）
subplot（2,1,2）
N=64;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）+0.8*randn（1,N）;y=dft（x,N）;
plot（q,abs（y））;title（'DFTN=64（withnoise）'）
由图可以看出这种噪音不影响信号检测。
4.对3，加大噪声到2*randn（1,N）和8*randn（1,N），画出并比较不同噪声下时域波形和频谱。
subplot（2,1,1）
N=64;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）+2*randn（1,N）;y=dft（x,N）;
plot（q,abs（y））;title（'DFTN=64（withnoise2）'）
subplot（2,1,2）
N=64;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）+8*randn（1,N）;y=dft（x,N）;
plot（q,abs（y））;title（'DFTN=64（withnoise8）'）
subplot（3,2,1）
N=64;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）+0.8*randn（1,N）;
plot（x）;title（'噪声为0.8*w的信号'）
y=dft（x,N）;
subplot（3,2,2）
plot（q,abs（y））;title（'噪声为0.8*w时的频谱'）
subplot（3,2,3）
N=64;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）+2*randn（1,N）;
plot（x）;title（'噪声为2*w时的信号'）
y=dft（x,N）;
subplot（3,2,4）
plot（q,abs（y））;title（'噪声为2*w时的频谱'）
subplot（3,2,5）
N=64;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）+8*randn（1,N）;
plot（x）;title（'噪声为8*w时的信号'）
y=dft（x,N）;
subplot（3,2,6）
plot（q,abs（y））;title（'噪声为8*w时的频谱'）
实验分析：
当噪声较小时，不影响信号的检测，但当噪声较大时，就看不出原信号的频率成分了，可以继续加大噪声，可看到其频谱杂乱无章了。
5.用一个N点DFT计算两个长度为N的实序列N点离散傅里叶变换，并将结果和直接使用两个N点DFT得到的结果进行比较。
x=[123456];
y=[654321];
[a,b]=dft_2（x,y）
a=
Columns1through3
21.0000-3.0000+5.1962i-3.0000+1.7321i
Columns4through6
-3.0000-3.0000-1.7321i-3.0000-5.1962i
b=
Columns1through3
21.00003.0000-5.1962i3.0000-1.7321i
Columns4through6
3.00003.0000+1.7321i3.0000+5.1962i
函数文件如下：
function[y1,y2]=dft_2（a,b）
N=length（a）;
x=zeros（1,N）;
x=a+j*b;
X=dft（x,N）;
X0=conj（fliplr（X））;
X0=[X0（N）X0（1:
N-1）];
y1=（X+X0）./2;
y2=（X-X0）./2./j;
直接运行计算：
dft（x）
ans=
Columns1through3
21.0000-3.0000+5.1962i-3.0000+1.7321i
Columns4through6
-3.0000-3.0000-1.7321i-3.0000-5.1962i
dft（y）
ans=
Columns1through3
21.00003.0000-5.1962i3.0000-1.7321i
Columns4through6
3.00003.0000+1.7321i3.0000+5.1962i
6．比较DFT和DFT的运算时间。
（计时函数tic，toc）
N分别取256，512，1024，2048，4096，…
程序如下：
N=256;
N=4096;
x=randn（1,N）;
tic
y=dft（x,N）;
toc
tic
z=dft（x）;
toc
N=256：
Elapsedtimeis0.172000seconds.
Elapsedtimeis0.015000seconds.
N=512：
Elapsedtimeis0.687000seconds.
Elapsedtimeis0.000000seconds.
N=1024：
Elapsedtimeis3.031000seconds.
Elapsedtimeis0.047000seconds.
N=2048：
Elapsedtimeis13.375000seconds.
Elapsedtimeis0.063000seconds.
N=4096：
Elapsedtimeis59.250000seconds.
Elapsedtimeis0.125000seconds.
7．对给定语音信号进行谱分析，写出采样频率，画出语音信号的波形及频谱，并分析语音信号的频率分布特点。
（1）画时域波形并对整个语音序列做DFT
[x,fs]=wavread（'C:
\ai1.wav'）;
subplot（2,1,1）
N=length（x）;
n=0:
N-1;
plot（n,x）;
xlabel（'n'）;
ylabel（'x'）;
title（'时域波形'）;
subplot（2,1,2）;
N=length（x）;
n=0:
N-1;
t=0.01*n;
q=n*2*pi/N;
y=dft（x,N）;
plot（q,abs（y））;
xlabel（'n'）;
ylabel（'ai1'）;
title（'DFT'）;
（2）并分别求出k=300，3500所对应的信号频率（Hz）
[x,fs]=wavread（'C:
\ai1.wav'）
N=length（x）;
n=0:
N-1;
t=n*（1/fs）;
q=n*2*pi/N;
n1=300;
q1=n1*2*pi/N;
f1=q1*fs/（2*pi）;
n2=3500;
q2=n2*2*pi/N;
f2=q2*fs/（2*pi）;
fs=
16000
（3）从语音中截取一段语音（256点）做DFT，得出频谱，画出时域波形及频谱。
分别求出k=5，60时所对应的信号频率（Hz）
程序如下：
[x,fs]=wavread（'C:
\ai1.wav'）;
subplot（2,1,1）;
N=256;
n=0:
N-1;
x=x（1:
256）;
plot（n,x）;
xlabel（'n'）;
ylabel（'x'）;
title（'256点时域波形'）;
subplot（2,1,2）;
N=256;
n=0:
N-1;
t=0.01*n;
q=n*2*pi/N;
x=x（1:
256）;
y=dft（x,N）;
plot（q,abs（y））;
xlabel（'n'）;
ylabel（'ai1'）;
title（'DFT-256'）;
（4）分别求出k=5，60时所对应的信号频率（Hz）的程序。
[x,fs]=wavread（'C:
\ai1.wav'）;
subplot（2,1,1）
N=256;
n=0:
N-1;
x=x（1:
256）;
t=0.01*n;
q=n*2*pi/N;
x=x（1:
256）;
y=dft（x,N）;
n1=5;
q1=n1*2*pi/N;
f1=q1*fs/（2*pi）;
n2=60;
q2=n2*2*pi/N;
f2=q2*fs/（2*pi）;
升级会员 