Matlab实现DFT

原文链接

前言：
摘自：
http://zhidao.baidu.com/link?url=lB3CS6qxlLIy7EuC8-X0MGNcFndC4ntrfXHdzdk5Wo_oJuRut8iHXoaOTEwA4BRBUeXpuibFuNd5Vp1hKsnsE_KrY8ALSESgHYpJHCPNKbO

　　MATLAB 傅里叶变换：
　　傅立叶变换的分类：
傅立叶级数：将周期性连续函数变换为离散频率点上的函数
傅立叶变换：将连续函数变换为连续频率的函数
离散时间傅立叶变换：将离散函数变换为连续频率的函数
离散傅立叶变换：将有限长离散函数变换为离散频率点上的函数
其中FFT是离散傅立叶变换的快速计算方法，适用于离散信号，并且注意变换后的点数与信号的采样点数一致。尽管可以将信号补0，但补0不能提高频域的分辨率。matlab中提供了函数fft做一维的FFT。

　　时域谱和频域谱是相互对应；时域的信号长度，决定频域的采样间隔，它们成导数关系；
　　时域中信号有N点，每点间隔dt，所以时域信号长度为N*dt；那么频谱每点的间隔F就是1/(N*dt)。（注：这里的采样间隔F=1/(N*dt)=Ws/N称为频谱分辨率，表示对x(n)在一个频谱周期内的N点等间隔采样）
　　傅立叶变换结果和原来信号有相同的点数，所以m=N，又第一点一定对应0频率，所以频域信号的很坐标就是(0:m-1)/(N*dt)，这句就是根据这个很坐标和频谱c，画出频谱plot((0:m-1)/(N*dt),c），所以在频谱图上，可以根据峰值的位置的横坐标读出对应的频率。
　　

clear all;
N=256;dt=0.02;
n=0:N-1;t=n*dt;
x=sin(2*pi*t);
m=N;
a=zeros(1,m);b=zeros(1,m);
for k=0:m-1
    for ii=0:N-1
        a(k+1)=a(k+1)+2/N*x(ii+1)*cos(2*pi*k*ii/N);
        b(k+1)=b(k+1)+2/N*x(ii+1)*sin(2*pi*k*ii/N);
    end
c(k+1)=sqrt(a(k+1)^2+b(k+1)^2);
end
subplot(211);plot(t,x);title('原始信号'),xlabel('时间/t');
f=(0:m-1)/(N*dt);
subplot(212);plot(f,c);hold on
title('Fourier');xlabel('频率/HZ');ylabel('振幅');
ind=find(c==max(c),1,'first');%寻找最大值的位置
%%%%%%% ind = find(X, k) 或ind = find(X, k, 'first')
%   返回第一个非零元素k的索引值（顺序）。
%   k必须是一个正数，但是它可以是任何数字数值类型。

x0=f(ind); %根据位置得到横坐标（频率）
y0=c(ind); %根据位置得到纵坐标（幅度）
plot(x0,y0,'ro');hold off

%%%%%%% hold on 和hold off，是相对使用的.通常是一个图上画两个曲线进行比较。
%   前者的意思是，你在当前图的轴（坐标系）中画了一幅图，再画另一幅图时，原来的图还在，与新图共存，都看得到
%   后者表达的是，你在当前图的轴（坐标系）中画了一幅图，此时，状态是hold off,则再画另一幅图时，原来的图就看不到了，在轴上绘制的是新图，原图被替换了

text(x0+1,y0-0.1,num2str(x0,'频率=%f'));

%   text(x,y,'string')在图形中指定的位置(x,y)上显示字符串string

   
   
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

---

栗子1：给x(n)绘图

global n;
n=-5:5;
x=sin(pi*n/5);
subplot(3,2,1);         %suplot定义fingure图形排列
%subplot（m，n，p）表示图排成m行，n列，p指的是你把曲线画到fingure中哪个图上，通常之后再加上plot(n,x);
stem(n,x,'.');          %将数据序列x从n值（即坐标（n，x）按照茎状形式画出并以'.'终止。
line([-5,6],[0,0]);     %line([起点横坐标，终点横坐标],[起点纵坐标，终点纵坐标]);
%表示画出从（-5,0）到（6,0）的一条直线，而不是（-5,6）到（0,0）.
axis([-5,6,-1.2,1.2]);
%axis([xmin xmax ymin ymax])来标注输出的图线的最大值与最小值，用来显示坐标范围。
%常与axis on/off来显示/关闭坐标轴上的标记，单位等
xlable('n');ylable('x(n)');
   
   
1
2
3
4
5
6
7
8
9
10
11
12

---

栗子2和栗子3分别对应书本P95和P98

---

栗子2:对待测信号xa(t)=(1+cos(2*pi*100*t)).*cos(2*pi*600*t)进行谱分析，若用fs=3kHz频率抽样，抽样点数为512

clear;
clc;
fs=3000;
N=512;
n=0:N-1;t=n/fs;
x=(1+cos(2*pi*100*t)).*cos(2*pi*600*t);

y=fft(x,N);
mag=abs(y);
f=n*fs/N;   %频率序列
subplot(2,1,1);plot(f,mag);
xlabel('频率/Hz');ylabel('振幅/X(k)');title('N=512,fs=3000HZ');grid on;

subplot(2,1,2);plot(1:N/2); %会出Nyquest频率之前随频率变化的振幅
xlabel('频率/Hz');ylabel('振幅/X(k)');title('Nyquest 图幅值');grid on;
%   grid on开启网格线

   
   
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

可知频率响应在f=500,600,700处突出，可断定该信号有三个f等于上的主要分量。

以下摘自：
http://wenku.baidu.com/link?url=x1hi17BxtAD_EskUAmK7wR5ppIypVFw7EW39YhN7jSZJnmP5o9Q485y9aJkl3brXX9wWIqWVRf-N5xy9MYcT4xJmi-CK0UhhYn6z22DeLW_

---

栗子3

%1.掌握DFT的概念和用法(DFT实质：有限长序列M傅里叶变换的有限点N的离散采样)
% 在频域分析信号分两种：
% （1）.对确定性信号进行傅里叶变换，分析频谱信息。
% （2）.随机信号的傅里叶信号不存在，转向研究它的功率谱。
% 功率谱是定义为时域信号傅氏变换模平方然后除以时间长度。是随机过程的统计平均概念，所表现的是单位频带内信号功率随频率的变换情况。
%2.利用DFT进行信号检测和谱分析

实验内容：

1．利用DFT计算信号功率谱。

t=0:0.001:6;
x=sin(2*pi*50*t)+sin(2*pi*120*t)+randn(1,length(t));
%   A = randn(m,n) 或 A = randn([m n])：返回一个m*n的随机项矩阵。
%   size(A)：获取数组A的行数和列数
%   length(A)：获取数组长度（即行数或列数中的较大值）
Y=fft(x,512);
P=Y.*conj(Y)/512;
%   '.*'运算符用于矩阵间对应元素的相乘，或数与数之间，数与矩阵之间的相乘。
%注：a*b表示矩阵a与矩阵b进行矩阵相乘。
%   a.*b表示矩阵a中的元素与矩阵b中的元素按照相同位置进行相乘，得到的结果作为新矩阵中相同位置的元素。
%   conj(A):返回A的共轭值
f=1000*(0:255)/512;
plot(f,P(1:256))
   
   
1
2
3
4
5
6
7
8
9
10
11
12
13

2.进行信号检测。分析信号频谱对应的数字频率和模拟频率之间的关系。模拟信号x(t)=2*sin(4*pi*t)+5*cos(8*pi*t),以t=0.01n（0<=n<=N-1）进行采样，求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=fft(x,N);
plot(q,abs(y));title('DFT N=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=fft(x,N);
plot(q,abs(y));title('DFT N=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=fft(x,N);
plot(q,abs(y));title('DFT N=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=fft(x,N);
plot(q,abs(y));title('DFT N=60')

%%%%%%%%%%%%注：X(k)为x(n)的DFT(x(ejw))在区间[0,2pi]上的N点等间隔采样
   
   
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

3.对2，进一步增加截取长度和DFT点数，如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=fft(x,N);
plot(q,abs(y));title('DFT N=256')
   
   
1
2
3
4
5

分析原因：采样频率f=100Hz，采样间隔T=0.01s，取样点数N数值越大时，则谱分辨率F=f/N越小，即分辨率越高。如当采样点数为N=256，即可区分两个频率不同的信号。

4.对3，加入噪声
取N=64，并在信号中加入噪声randn。

subplot(2,2,1);
f=100;
N=64;
n=0:N-1;t=n/f;
q=n*2*pi/N;
x=2*sin(4*pi*t)+5*cos(8*pi*t);
y=fft(x,N);
plot(q,abs(y));title('DFT N=64');

subplot(2,2,2);
x=2*sin(4*pi*t)+5*cos(8*pi*t)+0.8*randn(1,N);
y=fft(x,N);
plot(q,abs(y));title('DFT N=64(with 0.8 noise)');

subplot(2,2,3);
x=2*sin(4*pi*t)+5*cos(8*pi*t)+4*randn(1,N);
y=fft(x,N);
plot(q,abs(y));title('DFT N=64(with 4 noise)');

subplot(2,2,4);
x=2*sin(4*pi*t)+5*cos(8*pi*t)+16*randn(1,N);
y=fft(x,N);
plot(q,abs(y));title('DFT N=64(with 16 noise)');

   
   
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

可以看出，加大噪声并比较，可以看出噪声较小时，不影响信号检测，当噪声较大时，就看不出院信号的频率成分了。

---

加注：对上面进行DFT(横坐标是2*pi*n/N)和FFT(横坐标是fs*n/N)横坐标不同的解释。
摘自：
http://wenku.baidu.com/link?url=3DSg4R-sLHNJWIJh5n5QdNPaSzLOEA2VlnHcXYAj_NHFlzKvEP87sX9XA33CLQUtt4akNZrA9MHCcxXoqbemvUshcx4rMiw0LnX0V5ElKnm

DFT归一化的圆周频率为2*pi*n/N实际上是2*pi*F(F=f/fs),得后来的实际物理频率是fs*n/N。这也正是下文所介绍的从s域到z域的变换。
角频率（模拟频率）对应的是s平面，圆周频率对应的是z平面。
当模拟频率在s平面的虚轴上从0变到fs 时，数字频率在z平面单位圆上从0变到2 pi。这正是n/N然后一个乘以fs另一个乘以2*pi的原因。
当模拟频率在s平面的虚轴上从2fs变到4fs时，数字频率在z平面单位圆上仍然从0变到2 pi。因此后面的再研究就没有意义。