数字信号处理的MATLAB实现——快速傅里叶变换

一、FFT函数

1.调用格式:

X=fft(x,N)

(1)当x长度大于N时,以x序列的长度做FFT
(2)当x长度小于N时,在x后补零至N长再做FFT
(3)当没有设置N时,以x序列的长度做FFT

二、IFFT函数

1.调用格式

x=ifft(X)

三、计算初始相位角与初始幅度

1.理论基础

当f0与FFT后频谱上的某根谱线相重合时,可求出对应的信号赋值和初始相角。但如果f0在两条谱线之间,就不能用这个方法来计算,这就是栅栏效应。

2.调用格式

(1)计算初始相位角

theta=abs(X(k))

其中X(k)为第k点的FFT值
(2)计算初始幅度

A=2*abs(X(k)))/N

其中乘2是因为变换为单边谱,N为FFT点数,求出来的A为归一化幅度值。当频率在两条FFT谱线之间时,测量的幅度值有一定的误差。

3.范例


clear all; clc; close all;

fs=1000;                         % 采样频率
N=1000;                          % 信号长度
t=(0:N-1)/fs;                    % 设置时间序列
f1=50; f2=65.75;                 % 两信号频率
x=cos(2*pi*f1*t)+cos(2*pi*f2*t); % 设置信号
X=fft(x);                        % FFT
Y=abs(X)*2/1000;                 % 计算幅值
freq=fs*(0:N/2)/1000;            % 设置频率刻度
[A1, k1]=max(Y(45:65));          % 寻求第1个信号的幅值
k1=k1+44;                        % 修正索引号
[A2, k2]=max(Y(60:70));          % 寻求第1个信号的幅值
k2=k2+59;                        % 修正索引号
Theta1=angle(X(k1));
Theta2=angle(X(k2));
% 显示频率、幅值和初始相角
fprintf('f1=%5.2f   A1=%5.4f   Theta1=%5.4f\n',freq(k1),A1,Theta1); 
fprintf('f2=%5.2f   A2=%5.4f   Theta2=%5.4f\n',freq(k2),A2,Theta2);

% 作图
subplot 211; plot(freq,Y(1:N/2+1),'k'); xlim([0 150]); 
xlabel('频率/Hz'); ylabel('幅值'); title('频谱图');
subplot 223; stem(freq,Y(1:N/2+1),'k'); xlim([40 60]);
xlabel('频率/Hz'); ylabel('幅值'); title('50Hz分量');
subplot 224; stem(freq,Y(1:N/2+1),'k'); xlim([55 75]);
xlabel('频率/Hz'); ylabel('幅值'); title('65.75Hz分量');

四、频谱图中频率刻度的设置(以下讨论只考虑正频率)

1.N为奇数时

①freq=(0:(N-1)/2*fs/N)
③freq=linspace(0,fs/2-fs/N/2,(N-1)/2)

2.N为偶数时

①freq=(0:N/2-1)*fs/N
②freq=(0:N/2)*fs/N
③freq=linspace(0,fs/2,N/2+1)

3.范例

clear all; clc; close all;
fs=128;                         % 采样频率
N=128;                          % 信号长度
t=(0:N-1)/fs;                   % 时间序列
y=cos(2*pi*30*t);               % 余弦信号
y=fft(y,N);                     % FFT
%freq=(0:N/2-1)*fs/N;              % 按式(2-2-6c)设置正频率刻度
freq=linspace(0,fs/2,N/2)
% 作图
stem(freq,abs(y(1:N/2)),'k')
xlabel('频率(Hz)'); ylabel('幅值');
title('(b) 只有正频率刻度')
xlim([25 35]);
set(gcf,'color','w');

五、认识单频正弦信号相位谱

1.代码案例

%
% pr2_2_4 
clear all; clc; close all;
fs=1000;                       % 采样频率
f0=50;                         % 信号频率
A=1;                           % 信号幅值
theta0=pi/3;                   % 信号初始相角
N=1000;                        % 信号长度
t=(0:N-1)/fs;                  % 设置时间序列
x=A*cos(2*pi*f0*t+theta0);     % 设置信号
X=fft(x);                      % FFT
n2=1:N/2+1;                    % 设置索引号序列
freq=(n2-1)*fs/N;              % 设置频率刻度
% 第一部分
THETA=angle(X(n2));            % 计算初始相角
Am=abs(X(n2));                 % 计算幅值
ph0=THETA(51);                 % 计算信号的初始相角
fprintf('ph0=%5.4e   %5.4e   %5.4e\n',real(X(51)),imag(X(51)),ph0);
% 作图
%subplot 211; plot(freq,abs(X(n2))*2/N,'k');
%xlabel('频率/Hz'); ylabel('幅值')
%title('幅值谱图')
%subplot 212; plot(freq,THETA,'k')
%xlabel('频率/Hz'); ylabel('初始角/弧度')
%title('相位谱图')
%set(gcf,'color','w');
% 第二部分
Th=0.1;                        % 设置阈值
thetadex=find(Am<Th);          % 寻找小于阈值的那线谱线的索引
THETA1=THETA;                  % 初始化THETA1
THETA1(thetadex)=0;            % 对于小于阈值的那线谱线初始相位都为0
% 作图
figure(3)
%pos = get(gcf,'Position');
%set(gcf,'Position',[pos(1), pos(2)-100,pos(3),(pos(4)-160)]);
plot(freq,THETA1,'k')
xlabel('频率/Hz'); ylabel('初始角/弧度')
title('相位谱图')
set(gcf,'color','w');

2.分析

在这里插入图片描述
因为每个频率都有初始相位,所以相位谱非常乱。

3.解决方法

可以通过设置有用频率分量最小幅值的方法,筛选出无用频率分量,将其赋值置为0。得到的频谱图如下图所示:
在这里插入图片描述

六、消除趋势项

1.代码分析

clear all; clc; close all;
load qldata.mat                  % 读入数据
N=length(y);                     % 数据长度
time=(0:N-1)/fs;                 % 时间刻度
% 第一部分
Y=fft(y);                        % FFT
n2=1:N/2+1;                      % 取正频率索引序列
freq=(n2-1)*fs/N;                % 频率刻度
% 作图
subplot 211; plot(time,y,'k'); ylim([0 15]); grid;
title('有趋势项的数据')
xlabel('时间/s'); ylabel('幅值');
subplot 212; plot(freq,abs(Y(n2)),'k')
title('有趋势项的数据频谱')
xlabel('频率/Hz'); ylabel('幅值');
set(gcf,'color','w');
% 第二部分
x=detrend(y);                    % 消除趋势项
X=fft(x);                        % FFT
% 作图
figure
subplot 211; plot(time,x,'k'); ylim([-5 5]); grid;
title('消除趋势项后的数据')
xlabel('时间/s'); ylabel('幅值');
subplot 212; plot(freq,abs(X(n2)),'k');
title('消除趋势项后的数据频谱')
xlabel('频率/Hz'); ylabel('幅值');
set(gcf,'color','w');

2.问题分析

可以看到FFT频谱图中什么都没有,几乎都为0。这是因为波形中有一个明显的趋势项。正是因为该趋势项的存在使得频谱分析中有很大的直流分量。
在这里插入图片描述

3.解决方法

为了消除趋势项,可采用detrend函数:

xd=detrend(x)

消除趋势项后的FFT,如图所示:
在这里插入图片描述

七.将频谱图的纵坐标设置为分贝

1.方法:

(1)y轴用对数坐标

调用semilog函数

semilogx(x,y)  %x轴为对数坐标
semilogy(x,y)  %y轴为对数坐标
loglog(x,y)    %x,y轴都为对数坐标

(2)转换为分贝值后再作图

令Y=20*log10(abs(Y))

八、频谱的栅栏现象

在这里插入图片描述

九、频谱的泄露现象

1.定义

对于连续信号的采样序列进行DFT运算,由于时间长度取有限值,即将信号截断,使信号的带宽被扩展了,这种现象称为泄露。

2.通过加窗减小频谱泄漏

(1)如何选取窗函数

选择第一旁瓣衰减大,旁瓣峰值衰减快的窗函数有利于缓解截断过程中产生的频谱泄漏问题。

(2)N点DFT点数的选择

在这里插入图片描述
其中K为窗函数的主瓣宽度与矩形窗的主瓣宽度之比。

(3)理想窗函数

①窗函数的频谱尽可能窄,以提高谱估计时的频域分辨率和减小泄漏。
②尽量减小窗函数频谱的最大旁瓣的相对幅度,以使旁瓣高度随频率尽快衰减。

(4)解决加窗之后的幅值的衰减

为了修正加窗造成的幅值衰减,可在加窗FFT后的频谱进行修成,乘上窗函数的恢复系数
在这里插入图片描述

十、分辨率

1.定义

(1)物理分辨率:谱分析中将信号中两个靠的很近的谱峰仍能保持分辨的能力。
(2)计算分辨率使用DFT时,在频率轴上的所能得到的最小频率间隔。
在这里插入图片描述

十一、如何选择采样频率和信号长度

1.采样点数的选择

(1)示例:

在这里插入图片描述

(2)代码分析

①代码块:

clear all; clc; close all;
M=256; fs=10;                                   % 设置数据长度M和采样频率fs
f1=1; f2=2.5; f3=3;                             % 设置3个正弦信号的频率
t=(0:M-1)/fs;                                   % 设置时间序列
x=cos(2*pi*f1*t)+cos(2*pi*f2*t)+cos(2*pi*f3*t); % 计算出信号波形

X1=fft(x,20);                                   % FFT变换
X2=fft(x,40);
X3=fft(x,128);
freq1=(0:10)*fs/20;                              % 计算3个信号在频域的频率刻度
freq2=(0:20)*fs/40;
freq3=(0:64)*fs/128;
% 作图
plot(freq1,abs(X1(1:11)),'k',freq2,abs(X2(1:21)),'k',freq3,abs(X3(1:65)),'k');
%legend('N=20','N=40','N=128');
title('不同N值的DFT变换');
xlabel('频率/Hz'); ylabel('幅值');
set(gcf,'color','w');

②结果:N的点数不同,正弦信号的幅值也不同
在这里插入图片描述

十一.FFT中的补零问题

1.原理:

在FFT时经常会对输入序列补零,最常见的是补零后把N设置为2的k次幂。补零实际上是对FFT后的X(k)进行插值,补零后能改善栅栏效应,使原先不清晰的谱线显现出来。

2.代码分析

%
% pr2_2_12 
clear all; clc; close all;

fs=200;                            % 采样频率
f1=30; f2=65.5;                    % 两信号频率
N=200;                             % 信号长度
n=1:N;                             % 样点索引
t=(n-1)/fs;                        % 时间刻度
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);   % 信号

X1=fft(x);                         % 按N点进行FFT
freq1=(0:N/2)*fs/N;                % N点时正频率刻度
X1_abs=abs(X1(1:N/2+1))*2/N;       % 信号幅值

L=2*N;                             % 补零后FFT长度
X2=fft(x,L);                       % 按L长进行FFT
freq2=(0:L/2)*fs/L;                % L点时频率刻度
X2_abs=abs(X2(1:L/2+1))*2/N;       % 信号幅值
% 作图
subplot 211; plot(freq1,X1_abs,'k'); 
grid; ylim([0 1.2]);
xlabel('频率/Hz'); ylabel('幅值');
title('(a) 补零前FFT谱图')
subplot 212; plot(freq2,X2_abs,'k');
grid; ylim([0 1.2]);
xlabel('频率/Hz'); ylabel('幅值');
title('(b) 补零后FFT谱图')
set(gcf,'color','w');

由于在FFT上内插,所以可以看到主峰的周围有波纹震荡,这是由于泄露造成的。实际上补零前的信号并不是没有泄漏,而是在谱线的各个频点上正好与谱线的零点相重合。
在这里插入图片描述

在这里插入图片描述

3.注意:

补零不能增加物理分辨率,只能增加计算分辨率。要增加物理分辨率,唯一的方法是增加信号的有效长度。

十二、整周期采样

在这里插入图片描述

©️2020 CSDN 皮肤主题: 深蓝海洋 设计师:CSDN官方博客 返回首页