阅读量:0
文章目录
本文对matlab中fft的使用作出详细说明，并对频谱的双边、单边幅度谱与相位谱加以说明。
语法
Y = fft(X) Y = fft(X,N) Y = fft(X,N,dim) 说明
FFT是DFT的快速算法，当FFT点数为2的整数次幂时，MATLAB可以使用FFT的快速算法；如果不是2的整数次幂，那么只能使用公式的算法，实质上未使用上快速算法，两者在计算时间上有差异。
1.Y = fft(X) 用快速傅里叶变换 (FFT) 算法计算 X 的离散傅里叶变换 (DFT)。
- 如果
X是向量，则fft(X)返回该向量的傅里叶变换。即对于 Y = fft(X) 或 Y = fft(X,[],dim)，Y的大小等于 X 的大小。 - 如果
X是矩阵，则fft(X)将X的各列视为向量，并返回每列的傅里叶变换。 - 如果
X是一个多维数组，则fft(X)将沿大小不等于 1 的第一个数组维度的值视为向量，并返回每个向量的傅里叶变换。
2.Y = fft(X,N) 返回 N 点 DFT。即对于 Y = fft(X,n,dim)，size(Y,dim)的值等于 n，而所有其他维度的大小保持与在 X中相同。
- 如果
X是向量且X的长度小于N，则为X补上尾零以达到长度N。(下文有举例) - 如果
X是向量且X的长度大于N，则对X进行截断以达到长度N。(下文有举例) - 如果
X是矩阵，则每列的处理与在向量情况下相同。 - 如果
X为多维数组，则大小不等于 1 的第一个数组维度的处理与在向量情况下相同。
3.Y = fft(X,N,dim) 返回沿维度 dim 的傅里叶变换。例如，如果 X 是矩阵，则 fft(X,N,2) 返回每行的 N 点傅里叶变换。(下文有举例)
接下来对上述三种用法进行举例，实际绘制频谱中频谱分为单边谱和双边谱的绘制，在三种用法举例中顺带将这两种绘制方式一并列出了。
语法一：Y = fft(X)
fft(X)返回X长度的傅里叶变换
使用傅里叶变换求噪声中隐藏的信号的频率分量,其中使用单边谱的绘制方式呈现频谱图。
例：指定信号的参数，采样频率为 1 kHz，信号持续时间为 1.5 秒(即1500个信号点)。
clear all;close all;clc; %清理工作区，关闭所有窗口，清空文本 Fs = 1000; % Sampling frequency T = 1/Fs; % Sampling period L = 1500; % Length of signal t = (0:L-1)*T; % Time vector %构造一个信号，其中包含幅值为 0.7 的 50 Hz 正弦量和幅值为 1 的 150 Hz 正弦量。 S = 0.7*sin(2*pi*50*t) + sin(2*pi*150*t); X = S + 2*randn(size(t));%用均值为零、方差为 4 的白噪声扰乱该信号。 %% 绘制时域信号图 %在时域中绘制含噪信号。通过查看信号 X(t) 很难确定频率分量。 subplot(141); plot(t,X); title("Signal Corrupted with Zero-Mean Random Noise") xlabel("t (seconds)") ylabel("X(t)") subplot(143); plot(t,S); title("Original Signal") xlabel("t (seconds)") ylabel("X(t)") %% 计算加噪声信号的傅里叶变换 Y = fft(X);%此用法的fft点数为X的点数即1500 N = L;%fft点数为信号长度 %计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1。 P2 = abs(Y/N);%/N是对幅度的修正 P1 = P2(1:N/2+1); P1(2:end-1) = 2*P1(2:end-1);%第一个点为直流分量,是真实的幅值，保持不变；从第二个开始*2，由于DFT具有对称性看单边需要*2 f = (0:(N/2))*Fs/N;%其中Fs/N为频率分辨率 %定义频域 f 并绘制单侧幅值频谱 P1。与预期相符，由于增加了噪声，幅值并不精确等于 0.7 和 1。一般情况下，较长的信号会产生更好的频率逼近值。 subplot(142); plot(f,P1) title("Single-Sided Amplitude Spectrum of X(t)") xlabel("f (Hz)") ylabel("|P1(f)|") %% 现在，采用原始的、未破坏信号的傅里叶变换并检索精确幅值 0.7 和 1.0。 Y = fft(S); P2 = abs(Y/N); P1 = P2(1:N/2+1); P1(2:end-1) = 2*P1(2:end-1); subplot(144); plot(f,P1) title("Single-Sided Amplitude Spectrum of S(t)") xlabel("f (Hz)") ylabel("|P1(f)|") 使用fft(X)求得频谱绘制的单边谱如下：
语法二：Y = fft(X,N)
如果 X的长度小于 N，则为 X补上尾零以达到长度 N(FFT插值)
通过填充零来对信号的傅里叶变换进行插值。
例：指定信号的参数，采样频率为 80 Hz，信号持续时间为 0.8 秒。
clear all;close all;clc; %清理工作区，关闭所有窗口，清空文本 Fs = 80; T = 1/Fs; L = 65; t = (0:L-1)*T; %创建一个 2 Hz 正弦信号及其高次谐波的叠加。该信号包含一个 2 Hz 余弦波、一个 4 Hz 余弦波和一个 6 Hz 正弦波。 X = 3*cos(2*pi*2*t) + 2*cos(2*pi*4*t) + sin(2*pi*6*t); %在时域中绘制该信号。 subplot(131); figure(1); plot(t,X) title("Signal superposition in time domain") xlabel("t (ms)") ylabel("X(t)") %计算信号的傅里叶变换。 Y = fft(X);%不补零时 %计算信号的单侧幅值频谱。 f1 = Fs*(0:(L-1)/2)/L; P2 = abs(Y/L); P1 = P2(1:(L+1)/2); P1(2:end) = 2*P1(2:end); %在频域中绘制单侧频谱。由于信号的时间采样相当短，傅里叶变换的频率分辨率不够精确，不足以显示 4 Hz 附近的峰值频率。 subplot(132); plot(f1,P1,"-o") title("Single-Sided Spectrum of Original Signal") xlabel("f (Hz)") ylabel("|P1(f)|") %% 为了更好地评估峰值频率，您可以通过用零填充原始信号来增加分析窗的长度。这种方法以更精确的频率分辨率自动对信号的傅里叶变换进行插值。 %从原始信号长度确定是下一个 2 次幂的新输入长度。用尾随零填充信号 X 以扩展其长度。计算填零后的信号的傅里叶变换。 N = 2^nextpow2(L); Y = fft(X,N); %计算填零后的信号的单侧幅值频谱。由于信号长度 n 从 65 增加到 128，频率分辨率变为 Fs/n，即 0.625 Hz。 f2 = Fs*(0:(N/2))/N; P4 = abs(Y/L); P3 = P4(1:N/2+1); P3(2:end-1) = 2*P3(2:end-1); %绘制填零后的信号的单侧频谱。此新频谱在 0.625 Hz 的频率分辨率内显示 2 Hz、4 Hz 和 6 Hz 附近的峰值频率。 subplot(133); plot(f2,P3,"-o") title("Single-Sided Spectrum of Padded Signal") xlabel("f (Hz)") ylabel("|P3(f)|") 填充零来对信号的傅里叶变换进行插值效果对比如下图：
需要注意的是：
1.由于进行fft点数的不同，将双边谱转换为单边谱的计算有所区别。如代码中25行以及43行所示。
2.补零后，频率分辨率由原来的Fs/L变为Fs/N,补零以后能改善栅栏效应，使原先不清晰的谱线显现。虽然数据长度在补零后增长到N，但其有效长度还应该是L，且计算幅值是要以有效长度来计算的。参看代码23行和41行。参看FFT中的补零问题_fft补零的效果更加明显。
双边谱转换为单边谱
解释如下(注意MATLAB中的向量索引从1开始，以下解释为0开始)：
假设序列{y}的DFT序列{Y}长度也为N，表示为：
根据傅立叶变换的理论，这个序列中的后半部分实际上表示的是负频率信息。实序列的傅立叶变换的元素间存在共轭关系：
其中，常数s为：
这样，当N是奇数时，可以将序列{Y}表示为：
单边谱是长度为(N+1)/2的序列
当N是偶数时，序列{Y}可以表示为：
单边谱是长度为(N/2)+1的序列
此处权系数取法是w1=1，w2=2。
具体参看：信号处理：单边、双边频谱间的相互转换
作者拙见，供参考：DFT序列Y0为直流分量，余下的序列存在共轭关系。点数N为奇数时，除去Y0，余下序列为偶数均能成对共轭，各占了一半，转换为单边谱时乘2；点数N为偶数时，除去Y0，余下序列为奇数，中间空下一个不能成对被共用，不必乘2，成对的转换为单边谱需要乘2。
如果 X 的长度大于 N，则对 X 进行截断以达到长度 N。
如果 n 小于信号的长度，则 fft 忽略第 n 个条目之后的其余信号值，并返回截断后的结果。
使用傅里叶变换求两种频率分量拼的信号,其中使用双边谱的绘制方式呈现频谱图。
例：前1024个点为幅值为2，频率为50的信号；后1024个点为幅值为2，频率为100的信号；共2048个点。
fft(X,N)中的N取1000时可以看到频谱只有一个50的频率，由于是双边谱所以幅值各占一半。为什么只有一个频率啦，由于X信号的长度大于 N，所以将信号进行了截断，取前1000个点周期延拓求频谱。
clear all;close all;clc; %清理工作区，关闭所有窗口，清空文本 % 参数设置 fs = 1000; % 采样率 f_signal1 = 50; % 正弦信号频率1 f_signal2 = 100; % 正弦信号频率2 t1 = 0:1/fs:(1024-1)/fs; t2 = 1024/fs:1/fs:(2048-1)/fs; % 生成正弦信号 signal1 = 2*sin(2*pi*f_signal1*t1); signal2 = 2*sin(2*pi*f_signal2*t2); %合并两不同频率信号 t = [t1 t2]; signal = [signal1 signal2]; N = 2048; %fft点数 %计算fft频率轴,显示为频率f= n*(fs/N) f = (0:N-1)*fs/N; % 计算fftshift频率轴 fshift = (-N/2:N/2-1)*(fs/N);%注意范围 % 进行FFT fft_signal =(abs(fft(signal,N)))/N;%此用法设置FFT点数为N fft_signal_shift = fftshift(fft_signal); %% 绘制时域图和频谱图 figure(1); % 时域图 - 原始信号 subplot(1, 3, 1); plot(t, signal); title('Original Signal in Time Domain'); xlabel('Time (s)'); ylabel('Amplitude'); % 频谱图 - 原始信号 subplot(1, 3, 2); plot(f, fft_signal); title('Original Signal Spectrum(fft 2048)'); xlabel('Frequency (Hz)'); ylabel('Magnitude'); % 频谱图 - 原始信号 subplot(1, 3, 3); plot(fshift, fft_signal_shift); title('Original Signal Spectrum'); xlabel('Frequency (Hz)'); ylabel('Magnitude'); N取1000时的时域和频谱图：
N取2048时的频谱图：
可以看出由于频率分辨率和截断(此处相当于加矩形窗)带来的频谱泄露的问题。此时出现两个频率分量。
使用fft绘制的频谱图如下：
可以看出直接FFT的结果包含了各一半的两种频率分量。
语法三：Y = fft(X,N,dim)
dim — 沿其运算的维度 为正整数标量
沿其运算的维度，指定为正整数标量。如果不指定维度，则默认为第一个大于 1 的数组维度。
fft(X,[],1) 沿 X 的各列进行运算，并返回每列的傅里叶变换。
fft(X,[],2) 沿 X 的各行进行运算，并返回每行的傅里叶变换。
如果 dim 大于 ndims(X)，则 fft(X,[],dim) 返回 X。当指定 n 时，fft(X,n,dim) 将对 X 进行填充或截断，以使维度 dim 的长度为 n。
即dim=1，指定参数沿 X 的列使用 fft； dim =2 ，指定参数沿 X 的行使用 fft。(此处以dim =2 的行运算举例)：
clear all;close all;clc; %清理工作区，关闭所有窗口，清空文本 Fs = 1000; % Sampling frequency T = 1/Fs; % Sampling period L = 1000; % Length of signal t = (0:L-1)*T; % Time vector %创建一个矩阵，其中每一行代表一个频率经过缩放的余弦波。结果 X 为 3×1000 矩阵。第一行的波频为 50，第二行的波频为 150，第三行的波频为 300。 x1 = cos(2*pi*50*t); % First row wave x2 = cos(2*pi*150*t); % Second row wave x3 = cos(2*pi*300*t); % Third row wave X = [x1; x2; x3];%X为3*1000double型 %在单个图窗中按顺序绘制 X 的每行的前 100 个项，并比较其频率。 figure(1); for i = 1:3 subplot(3,1,i) plot(t(1:100),X(i,1:100)) title("Row " + num2str(i) + " in the Time Domain") end %指定 dim 参数沿 X 的行（即对每个信号）使用 fft。 dim = 2; %计算信号的傅里叶变换。 Y = fft(X,L,dim); %计算每个信号的双侧频谱和单侧频谱。 P2 = abs(Y/L); P1 = P2(:,1:L/2+1); P1(:,2:end-1) = 2*P1(:,2:end-1); %在频域内，为单个图窗中的每一行绘制单侧幅值频谱。 figure(2); for i=1:3 subplot(3,1,i) plot(0:(Fs/L):(Fs/2-Fs/L),P1(i,1:L/2)) title("Row " + num2str(i) + " in the Frequency Domain") end 时域波形：
单边幅值频谱图：
相位
创建一个由频率为 15 Hz 和 40 Hz 的两个正弦波组成的信号。第一个正弦波是相位为 −π/4 的余弦波，第二个正弦波是相位为 π/2 的余弦波。以 100 Hz 的频率对信号进行 1 秒钟的采样。
clear all;close all;clc; %清理工作区，关闭所有窗口，清空文本 Fs = 100; t = 0:1/Fs:1-1/Fs; x = cos(2*pi*15*t - pi/4) + cos(2*pi*40*t + pi/2); %计算信号的傅里叶变换。将变换幅值绘制为频率函数。 y = fft(x); z = fftshift(y); figure(1); ly = length(y); f = (-ly/2:ly/2-1)/ly*Fs; stem(f,abs(z)) title("Double-Sided Amplitude Spectrum of x(t)") xlabel("Frequency (Hz)") ylabel("|y|") grid %计算变换的相位，删除小幅值变换值。将相位绘制为频率函数。 tol = 1e-6; z(abs(z) < tol) = 0;%删除小幅值变换值 theta = angle(z); figure(2); stem(f,theta/pi) title("Phase Spectrum of x(t)") xlabel("Frequency (Hz)") ylabel("Phase/\pi") grid 幅度和相位谱：
补充
向量的离散傅里叶变换
Y = fft(X) 和 X = ifft(Y) 分别实现傅里叶变换和傅里叶逆变换。对于长度为 n 的 X 和 Y，这些变换定义如下：
频谱泄露的解释
文章FFT中的补零问题_fft补零3.1 频谱泄漏仿真分析处及之后部分。
参考：
