欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识！
您现在的位置是： 首页

FFT：信号的频谱分析

程序员文章站 2022-03-10 16:22:55
...

对于下面这句话该怎么理解？

假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。

思考：一个序列波形，怎样求得峰值及其位置，进而求得周期？

（注：要先进行平滑滤波处理。）

摘自：http://www.ilovematlab.cn/thread-289950-1-1.html

前言：
1）可用[a,b]=max(data)(data是一离散数列)求极大值，但这种方法作用有限，不能找到连续极值。
2）[a,b] = findpeaks(data)：a 对应峰值，b 对应峰值位数。但缺点是只能找波峰值，无法找波谷值。

1.自建一个求峰值坐标函数

function A=zhao(M);
global n;
double M;double A;
n=1;
l=length(M);
a=zeros(1,(l-1));b=1:l-1;
A=[b' a' a'];       %对A进行初始化，(n-1)*3型
for i=1:l-1;
    if(M(i)>M(i+1)&&M(i-1))
        A(n,2)=i;A(n,3)=M(i);
        n=n+1;
    end
end
end
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14

FFT：信号的频谱分析

继而再求出周期（对于循环序列波形来说）：也就是相邻两个极大值之间的距离。

只需要再提取B的第二列（最高点出现的横坐标位置），做差分。若有负值，将负值及负值之后的全0数据都舍去。然后再将负值前面的s个数据相加再/s即可。

m=(B(:,2));     
p=diff(m);
t=0;
for i=1:l-2
    if(p(i)>0)
        t=p(i)+t;
        s=i;
    end

T=t/s;
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10

2.添注：用matlab函数找序列波形的零点

前言：基本波形绘制

%%%%%%%%绘制在-2pi~2pi区间冲激函数y=sin(t)/t的时域和频域图
clear;
n=400;
delta=4*pi/n;
t=-2*pi:delta:2*pi;
y=sinc(t);      %定义冲激函数
subplot(121);
plot(t,y);axis([-7,7,-0.4,1.3]);
title('sinc(t)');
grid on;
subplot(122);
Y1=fft(y,2048);
Y=fftshift(Y1);
c=-1024:1023;
plot(c,abs(Y));
axis([-500,500,-5,40]);
title('sinc(t)频谱');
grid on;
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18

FFT：信号的频谱分析

对于fftshift（）；

解释：摘自http://blog.sina.com.cn/s/blog_68f3a4510100qvp1.html

第三段：根据Fourier分析的相关结论，我们知道时域的采样将会造成频域的周期化，该周期为采样频率f_s（即书本上的Fs），我们谱分析是就是观察0~Fs这一段。

由奈奎斯特采样定理可以知道，fs必须≥信号最高频率的2倍才不会发生信号混叠，因此fs能采样到的信号最高频率为fs/2.所以只有f=fs/2范围内的信号才是被采样到的有效信号。那么，在w的范围内，得到的频谱肯定是关于n/2对称的。（所以下文的频谱仪显示窗的频率范围为0~Fs/2，是为了便于观察和分析）

和后面：如果我们不使用fftshift，其变换后的横坐标为0:f_s/(N-1):f_s， 如果使用fftshift命令，0频率分量将会移到坐标中心，这也正是matlab中帮助中心给出的意思：对fft的坐标进行了处理。实际上由于频谱的周期性，我们这样做是合理的，可以接受的。

一.频谱仪的时窗谱线分析

频谱仪用FFT完成数据流从时域到频域的变换。FFT长度为N(一般为2的n次幂)，N的大小，即时窗的长短，决定着频谱仪的频率分辨率，时窗越长，分辨率越高，可以将相隔很近的谱线区分开。

1.用矩阵运算的方法计算方波的DFS（离散傅里叶级数），在60的时窗宽度上方波宽度分别为5和12，并且画出x（n）和DFS（x（n））的杆状图。

clc;
L=5;N=60;
k=[-N/2:N/2];
xn=[zeros(1,(N-L+1)/2),ones(1,L),zeros(1,(N-L-1)/2)];    %方波信号的另一种表示
n=[0:N-1];

subplot(221);stem(n,xn);grid;
axis([0,60,-0.3,1.3]);title('xn');

p=[0:N-1];
WN=exp(-j*2*pi/N);
nk=n'*p;
WNnk=WN.^nk;
Xk=xn*WNnk;
magXk=abs([Xk(N/2+1:N) Xk(1:N/2+1)]);

%   这一段是求和，好好理解（我也不会....）

subplot(222);stem(k,magXk);grid
axis([-N/2,N/2,-0.5,5.5]);
xlabel('k');ylabel('spectrum');
title('DFS:L=5,N=60');grid;
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22

FFT：信号的频谱分析

注：方波的宽度越窄，响应频谱的带宽越宽。如冲激信号的频谱在无线区间上是一恒定值

2.+不同窗函数的频谱分析

a.

%%%%不同时窗正弦信号频谱的分析
%   设一个1Hz的正弦波，取10和50个周期，都做4096点FFT，比较二者的频谱
clear all;
t=0:0.1:10;
y=sin(2*pi*t);
Y1=fft(y,4096);
Y=fftshift(Y1);
c=[0:2047]/409.6;
subplot(121);
plot(c,abs(Y(2049:4096)));      %谁能告诉我这三行是什么意思吗？如为什么要/409.6，还有y的取值
axis([0,2,-5,60]);
title('sin(t)周期T=10');
grid;

t=0:0.1:50;
y=sin(2*pi*t);
Y1=fft(y,4096);
Y=fftshift(Y1);
c=[0:2047]/409;
subplot(122);
plot(c,abs(Y(2049:4096)));
axis([0,2,-5,60]);
title('sin(t)周期T=50');
grid;
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24

FFT：信号的频谱分析

解释：无限长的正弦信号频谱是一根根垂直于频率轴的垂直线段，有限长的正弦信号等于无限长正弦信号与时窗等长且高度为1的方波的乘积，由卷积定理可知，有限长的正弦信号的频谱等于1/(2*pi)方波的频谱与正弦信号的谱线的卷积。而窄方波的频谱宽，宽方波的频谱窄，故将t取值为50可得到更宽时窗的频谱，且失真较小。

b.加窗函数后的频谱分析。


%%%   设一个1Hz的正弦波，取20个周期，求加上窗函数（汉明窗）之后的频谱

clear all;
t=0.1:0.1:20;
y=sin(2*pi*t);
w=hamming(200);
y1=y.* w';
y2=fft(y,4096);
y3=fftshift(y2);
c=[0:2047]./409.6;
plot(c,abs(y3(2049:4096)));
axis([0,2,-5,100]);
title('sin(t)加汉明窗后频谱');
grid on;
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15

FFT：信号的频谱分析

解释：求频谱加不加函窗数的区别：不加窗函数是正弦谱线和矩形窗频谱的卷积，加窗函数是正弦谱线和汉明窗频谱的卷积。

二.周期序列的DFS及DFT

a.基本方波序列

FFT：信号的频谱分析

%%%%求序列x(n)=RN(n)的X(e.^jw)，并画出序列、模、相角曲线
%矩形序列：起点n0，终点n1，序列为1值区间[ns,nf].

clear all;
clc;
n0=input('n0=');n1=input('n1=');ns=input('ns=');nf=input('nf=');
if(n0>=n1|ns>=nf|n0>=ns|nf>=n1)
    error('输入错误');
end;
n=n0:n1;
x=[zeros(1,ns-n0),ones(1,nf-ns+1),zeros(1,n1-nf)];      %单位阶跃序列产生
subplot(311);stem(n,x);title('矩形序列RN(n)');
axis([-5,10,-0.2,1.2]);grid;

k=-500:1200;w=(pi/500)*k;       %[0,pi]段分为501点
X=x*(exp(-j*pi/500)).^(n'*k);       % 用矩阵-向量乘法求DTFT   

subplot(312);plot(w,abs(X));grid;
xlabel('');ylabel('模值');title('模值部分');
subplot(313);plot(w,angle(X));grid;
xlabel('以pi为单位的频率');ylabel('弧度');title('相角部分');
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21

FFT：信号的频谱分析
FFT：信号的频谱分析

b.周期序列

FFT：信号的频谱分析

条件同上，将其以N=8为周期进行周期延拓，得到序列x(n)~.画出DTFT.

%%%%求序列x(n)=RN(n)的X(e.^jw)，并画出序列、模、相角曲线
%矩形序列：起点n0，终点n1，序列为1值区间[ns,nf].

clear all;
clc;
n0=input('n0=');n1=input('n1=');ns=input('ns=');nf=input('nf=');
if(n0>=n1|ns>=nf|n0>=ns|nf>=n1)
    error('输入错误');
end;
n=n0:n1;
x=[zeros(1,ns-n0),ones(1,nf-ns+1),zeros(1,n1-nf)];      %单位阶跃序列产生

xtidle=x'*ones(1,n1-n0+1);  %用x转置乘以x的长度大小的全1矩阵。
xtidle=xtidle(:);          % A(:)就是按matlab中的存储顺序，从A(1)到A(end)依次按列排序,第二列接到第一列上，然后依次反复，使A变成一个列向量
xtidle=xtidle';     再转置形成1*（n1-n0+1）^2行矩阵
%%%%%%%%%%%%上面三行就是进行周期延拓的

n=0:length(xtidle)-1;

subplot(311);stem(n,xtidle);title('矩形序列RN(n)');
axis([0,80,0,1.1]);grid;

k=-500:1200;w=(pi/500)*k;       %[0,pi]段分为501点
X=xtidle*(exp(-j*pi/500)).^(n'*k);       % 用矩阵-向量乘法求DTFT   

subplot(312);plot(w,abs(X));grid;
xlabel('');ylabel('模值');title('模值部分');
subplot(313);plot(w,angle(X));grid;
xlabel('以pi为单位的频率');ylabel('弧度');title('相角部分');
  • 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

FFT：信号的频谱分析

三.FFT实现周期信号的频谱分析（暂略，待补）

上一篇： 02-方法引用、构造器引用、数组引用

下一篇： imx6ull移植mplayer

推荐阅读