再续FFT：信号的频谱分析

摘自：http://www.ilovematlab.cn/thread-119939-1-1.html
http://www.360doc.com/content/13/1208/18/13670635_335496776.shtml

---

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

假设采样频率为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

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

只需要再提取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;

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;

---

对于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;

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

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的方波的乘积，由卷积定理可知，有限长的正弦信号的频谱等于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;

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

二.周期序列的DFS及DFT

a.基本方波序列

%%%%求序列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('相角部分');

b.周期序列

条件同上，将其以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('相角部分');

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