声信号处理工程研究（2）---频域图（频谱）

引言：

开始之前先了解一下声信号的三元素（我自认为确定这三元素，就唯一确定了一个单频信号）。

针对一个单频信号，确定了幅值A，频率f，初相位φ（大概意思就是采样数据的第一个采样点对应的三角函数的角度）后，可以构建单频信号的模拟信号方程,也就唯一确定了一个单频信号，

如y（t）=A*cos(2*π*f*t+φ)。

解释一下相位：(2πft+φ)的值，称做 相位<==>相角，  φ称为初相<==>初相位<==>初相角。

初相为0时，100Hz的单频幅值为1的余弦波的波形图如下，刚好是一个简单的余弦函数图

初相为π/2时，100Hz的单频幅值为1的余弦波的波形图如下，变成正弦函数图了，0点的值应该为0，这里为3.062e-17是软件计算误差的问题，其他也是无限接近0

注：一般声信号使用cos函数，sin函数和cos函数其实一样，只是相位差π/2

正文：

什么是频域？

频域：自变量是频率,即横轴是频率,纵轴是该频率信号的幅度（声信号中是声压）,也就是通常说的频域图。

什么是频域图？

频域图：描述了信号的频率结构 及 频率与该频率信号幅度的关系。

处理流程是啥？

一个模拟信号（即现实中的声音），经过ADC采样之后，就变成了数字信号（某个音频文件，或采样数据数组等）。根据采样定理，采样频率要大于信号频率的两倍。

采样得到的数字信号，就可以做FFT变换了。N个采样点，经过FFT之后，就可以得到N个点的FFT结果（复数数组，每个复数为一个对应频率的幅值信息）。为了方便进行FFT运算，通常N取2的整数次方.

注：MATLAB调用更高级的快速傅里叶变换库，可以对任意基数(不一定非要2的整数次方)进行FFT,如100个采样点输入，经过FFT后，可以得到100个点的FFT结果。

假设采样频率为fs，信号频率f，采样点数为N。那么FFT之后结果就是一个为N个点的复数数组。每一个点就对应着一个频率点的信息。这个点的模值，就是该频率值下的幅值特性。

这个点的复数与实数的正切函数值，就是在该频率下的信号的相位。

怎么得到原始信号的幅值？

假设原始单频信号为A*cos(2*π*f*t+φ)，则幅值为A，那么FFT的得到的数组的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。

而第一个点就是直流分量（即0频率，零频分量），它的模值就是直流分量的N倍。

第一个点表示直流分量（即0Hz），而最后一个点N的再下一个点（实际上这个点是不存在的，这里是假设的第N+1个点，也可以看做是将第一个点分做两半分，另一半移到最后）则表示采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。

例如某点n所表示的频率为：Fn=(n-1)*Fs/N，即FFT得到数组X[N]，X[1]的模值为直流分量的信息，X[2]为(2-1)*Fs/N频率的幅值信息，N=1,2...N;

Fn所能分辨到频率为为Fs/N，Fs/N也叫频率分辨率。

问题/坑：频率分辨率越小，是不是能分辨更细，更多的频率的幅值信息？如果是，怎么减小Fs/N？

如果采样频率Fs为1024Hz，采样点数N为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到0.5Hz。如果要提高频率分辨力，则必须增加采样点数，也即采样时间。频率分辨率和采样时间是倒数关系。实际工程应用中不现实，需要在较短的时间内完成分析。

假设FFT之后某点n用复数a+bi表示，那么这个复数的模就是An=sqrt(a*a+b*b)，相位就是Pn=atan2(b,a)。根据以上的结果，就可以计算出n点（n≠1，且n<=N/2）对应的信号的表达式为：An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。对于n=1点的信号，是直流分量，幅度即为A1/N。
由于FFT结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。(后续会解释)

实战

下面以一个实际的信号来做说明。

针对单频模拟信号，初相为0，100Hz的单频幅值为1.4的余弦波

clc,clear;close all;
f=100;%设定模拟信号的频率为100Hz
fs=10000;%采样频率为10000Hz
t=0:1/fs:9999/fs;%生成采样时刻数组t      
data=1.4*cos(2*pi*f*t);%得到采样数据数组data，构建原始模拟信号 幅值根号2 ，频率为100hz ，初始相位为 0度 ，
%画时域图                  
figure(1);%创建图窗1
plot(t,data);
title('时域图');
xlabel('时间/s');
ylabel('振幅');
xlim([0 0.1]);

%画频域图 
figure(2);%创建图窗2
n=length(data);%采样数据长度n;
f=0:fs/n:fs/n*(n-1);%频率数组，频域图x轴
y=fft(data);% y为频域双边频域复数数组，不是频域图幅值数组；
        % 默认fft计算点数为采样数据data长度，fft计算点数=输出数组y的长度
y1=abs(y);%取模得到，模值数组
y21=y1/n*2;%模值除以n再乘2,得到频域图单边频谱幅值数组,也可写成y21=y1/（n/2）
y21(1)=y1(1)/2;%第一个值为0频分量，即为直流分量的幅值，真实幅值应为 y1(1)/n，这里将它修正为真实值
plot(f(1:n/2+1),y21(1:n/2+1));%画出频域图单边谱，fft只能得到采样率一半的频率的信息
title('频域图');
xlabel('频率/Hz');
ylabel('振幅');


%画相位图 
figure(3);%创建图窗3
p=angle(y);%相位数组，相位图y轴,弧度制
p1=angle(y)*180/pi;%相位数组，相位图y轴，角度制
plot(f(1:n/2+1),p1(1:n/2+1));
title('相位图');
xlabel('频率/Hz');
ylabel('相位/度');

时域图为

频域图为

可见，100Hz处有明显的谱线或者峰值，为模拟信号的频率与幅值，表示针对采集信号进析FFT后，得到了采集信号的频率信息，，推断采集信号是一个单频信号，其他频率存在值，但是和100HZ幅值差别很大（fft计算存在误差的原因）

相位图为

理论上应该只有100Hz有对应角度，没有其他频率的声音，所以其他频率相位为0，但是上面说了fft计算存在误差，其他频率也会计算出值，但是很小，复数的模接近于0，所以计算其他频率对应fft的复数值也会计算出角度。

针对混合信号，比如两个频率的混合声音，将初相为0，100Hz的单频幅值为1.4的余弦波和初相为0，500Hz的单频幅值为1.4的余弦波混合

clc,clear;close all;
f1=100;%设定模拟信号1的频率为100Hz
f2=500;%设定模拟信号2的频率为500Hz
fs=10000;%采样频率为10000Hz
t=0:1/fs:9999/fs;%生成采样时刻数组t      
data=1.4*cos(2*pi*f1*t)+1.4*cos(2*pi*f2*t);%得到采样数据数组data，构建原始模拟信号
%画时域图                  
figure(1);%创建图窗1
plot(t,data);
title('时域图');
xlabel('时间/s');
ylabel('振幅');
xlim([0 0.1]);

%画频域图 
figure(2);%创建图窗2
n=length(data);%采样数据长度n;
f=0:fs/n:fs/n*(n-1);%频率数组，频域图x轴
y=fft(data);% y为频域双边频域复数数组       
y1=abs(y);%取模得到，模值数组
y21=y1/n*2;%模值除以n再乘2,得到频域图单边频谱幅值数组,也可写成y21=y1/（n/2）
y21(1)=y1(1)/2;%第一个值为0频分量，即为直流分量的幅值
plot(f(1:n/2+1),y21(1:n/2+1));%画出频域图单边谱，fft只能得到采样率一半的频率的信息
title('频域图');
xlabel('频率/Hz');
ylabel('振幅');

时域图

频域图

思考：FFT的缺陷或者频域图的局限性？

比如：第一秒只有100Hz信号，后一秒只有500Hz信号的采集信号与第一秒只有500Hz信号，后一秒只有100Hz信号的采集信号两种不同信号分析得到的频域图是相同的。

%%%%%%%%%%%%%%%FFT与STFT%%%%%%%%%%%%%%%%%%%%%%%%%%
f1=100;%设定模拟信号1的频率为100Hz
f2=500;%设定模拟信号2的频率为500Hz
fs=10000;%采样频率为10000Hz

t1=0:1/fs:999/fs;     
data1=1.4*cos(2*pi*f1*t1);%得到采样数据数组data，构建原始模拟信号
data2=1.4*cos(2*pi*f2*t1);
data11=[data1 data2];%拼接两个数组，构建模拟信号，表示前一秒只有100Hz信号，后一秒只有500Hz信号
data22=[data2 data1];%拼接两个数组，构建模拟信号，表示前一秒只有500Hz信号，后一秒只有100Hz信号
t=0:1/fs:1999/fs;     %采样时刻数组t 
n=length(data11);%采样数据长度n;
%画时域图                  
figure(1);%创建图窗1
subplot(2,2,1),plot(t,data11);title('时域图1');xlabel('时间/s');ylabel('振幅');%xlim([0 0.1]);
subplot(2,2,3),plot(t,data22);title('时域图2');xlabel('时间/s');ylabel('振幅');%xlim([0 0.1]);
%画频域图 


f=0:fs/n:fs/n*(n-1);%频率数组，频域图x轴
y11=fft(data11);% y为频域双边频域复数数组，不是频域图幅值数组；
        % 默认fft计算点数为采样数据data长度，fft计算点数=输出数组y的长度
y1=abs(y11);%取模得到，模值数组
y21=y1/n*2;%模值除以n再乘2,得到频域图单边频谱幅值数组,也可写成y21=y1/（n/2）
y21(1)=y1(1)/2;%第一个值为0频分量，即为直流分量的幅值，真实幅值应为 y1(1)/n，这里将它修正为真实值
subplot(2,2,2),plot(f(1:n/2+1),y21(1:n/2+1));title('频域图1');xlabel('频率/Hz');ylabel('振幅');

y22=fft(data22);     
y2=abs(y22);%取模得到，模值数组
y22=y1/n*2;%模值除以n再乘2,得到频域图单边频谱幅值数组,也可写成y21=y1/（n/2）
y22(1)=y1(1)/2;%第一个值为0频分量，即为直流分量的幅值，真实幅值应为 y1(1)/n，这里将它修正为真实值
subplot(2,2,4),plot(f(1:n/2+1),y22(1:n/2+1));title('频域图2');xlabel('频率/Hz');ylabel('振幅');

参考链接：FFT的详细解释，相信你看了就明白了