前言：

摘自：

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)上显示字符串string1

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。因此后面的再研究就没有意义。