matlab - MATLAB fft 与数学傅立叶
问题描述
我正在使用以下代码生成信号的 fft 和数学傅立叶变换。然后我想在数学上重新创建 fft 的原始信号。这适用于数学信号,但不适用于 fft,因为它是离散变换。有谁知道我可以对我的逆变换方程进行哪些更改以使其适用于 fft?
clear all; clc;
N = 1024;
N2 = 1023;
SNR = -10;
fs = 1024;
Ts = 1/fs;
t = (0:(N-1))*Ts;
x = 0.5*sawtooth(2*2*pi*t);
x1 = fft(x);
Magnitude1 = abs(x1);
Phase1 = angle(x1)*360/(2*pi);
for m = 1:1024
f(m) = m; % Sinusoidal frequencies
a = (2/N)*sum(x.*cos(2*pi*f(m)*t)); % Cosine coeff.
b = (2/N)*sum(x.*sin(2*pi*f(m)*t)); % Sine coeff
Magnitude(m) = sqrt(a^2 + b^2); % Magnitude spectrum
Phase(m) = -atan2(b,a); % Phase spectrum
end
subplot(2,1,1);
plot(f,Magnitude1./512); % Plot magnitude spectrum
......Labels and title.......
subplot(2,1,2);
plot(f,Magnitude,'k'); % Plot phase spectrum
ylabel('Phase (deg)','FontSize',14);
pause();
x2 = zeros(1,1024); % Waveform vector
for m = 1:24
f(m) = m; % Sinusoidal frequencies
x2 = (1/m)*(x2 + Magnitude1(m)*cos(2*pi*f(m)*t + Phase1(m)));
end
x3 = zeros(1,1024); % Waveform vector
for m = 1:24
f(m) = m; % Sinusoidal frequencies
x3 = (x3 + Magnitude(m)*cos(2*pi*f(m)*t + Phase(m)));
end
plot(t,x,'--k'); hold on;
plot(t,x2,'k');
plot(t,x3,'b');```
解决方案
关于傅里叶变换有一些评论,我希望我能为你解释一切。另外,我不知道您所说的“数学傅立叶变换”是什么意思,因为您的代码中没有一个表达式类似于锯齿波的傅立叶级数。
要准确理解fft 函数的作用,我们可以一步一步做。首先,按照您的代码,我们创建并绘制一个周期的锯齿波。
n = 1024;
fs = 1024;
dt = 1/fs;
t = (0:(n-1))*dt;
x = 0.5*sawtooth(2*pi*t);
figure; plot(t,x); xlabel('t [s]'); ylabel('x');
我们现在可以计算一些东西。首先,奈奎斯特频率,即样本中可检测到的最大频率。
f_max = 0.5*fs
f_max = 512
此外,最小可检测频率,
f_min = 1/t(end)
f_min = 1.000977517106549
现在用 MATLAB 函数计算离散傅里叶变换:
X = fft(x)/n;
该函数获得离散傅里叶变换每一项的复系数。请注意,它使用 exp 表示法计算系数,而不是根据正弦和余弦。除以n
是为了保证第一个系数等于样本的算术平均值
如果要绘制变换信号的幅度/相位,可以键入:
f = linspace(f_min,f_max,n/2); % frequency vector
a0 = X(1); % constant amplitude
X(1)=[]; % we don't have to plot the first component, as it is the constant amplitude term
XP = X(1:n/2); % we get only the first half of the array, as the second half is the reflection along the y-axis
figure
subplot(2,1,1)
plot(f,abs(XP)); ylabel('Amplitude');
subplot(2,1,2)
plot(f,angle(XP)); ylabel('Phase');
xlabel('Frequency [Hz]')
这个情节是什么意思?它在图中显示了代表原始信号(锯齿波)的傅里叶级数项的复系数的幅度和相位。您可以使用此系数根据(截断的)傅里叶级数获得信号近似值。当然,要做到这一点,我们需要整个变换(不仅仅是前半部分,因为通常会绘制它)。
X = fft(x)/n;
amplitude = abs(X);
phase = angle(X);
f = fs*[(0:(n/2)-1)/n (-n/2:-1)/n]; % frequency vector with all components
% we calculate the value of x for each time step
for j=1:n
x_approx(j) = 0;
for k=1:n % summation done using a for
x_approx(j) = x_approx(j)+X(k)*exp(2*pi*1i/n*(j-1)*(k-1));
end
x_approx(j) = x_approx(j);
end
注意:上面的代码是为了澄清,并不打算很好地编码。求和可以在 MATLAB 中以比使用 for 循环更好的方式完成,并且代码中会弹出一些警告,警告用户为速度预先分配每个变量。
上面的代码使用截断的傅立叶级数计算x(ti)
每次的 。ti
如果我们同时绘制原始信号和近似信号,我们得到:
figure
plot(t,x,t,x_approx)
legend('original signal','signal from fft','location','best')
原始信号和近似信号几乎相等。事实上,
norm(x-x_approx)
ans = 1.997566360514140e-12
几乎为零,但不完全为零。此外,由于在计算近似信号时使用了复系数,上图将发出警告:
警告:复杂 X 和/或 Y 参数的虚部被忽略
但是您可以检查虚数项是否非常接近于零。由于计算中的舍入误差,它并不完全为零。
norm(imag(x_approx))
ans = 1.402648396024229e-12
请注意在上面的代码中如何解释和使用 fft 函数的结果以及它们如何以 exp 形式表示,而不是像您编码的那样根据正弦和余弦表示。
推荐阅读
- substrate - 升级到 Frame v2 版本
- c# - 来自父项控件的 WPF 对象绑定
- delphi - 通过 skSince/skSentSince 搜索 IMAP 邮箱会忽略 TDateTime 的时间部分
- python - 我做了一个程序,它通过输入 IP 地址来告诉位置等,但我现在必须这样格式显示错误,我应该怎么做才能纠正它
- python - 视图函数 Flask 中的额外参数
- android - 在我旋转屏幕之前,LiveData 不会加载到片段中
- javascript - 收到并处理完事件后需要关闭MongoDB Changestream吗?
- ibm-mq - 我无法将我的主要安装更改为 opt/mqm/
- react-native - 反应原生步进指示器
- rust - 是否可以忽略 Rust 的指数?