matlab自带的fft
没用明白,干脆自己写一个,牺牲了时间复杂度,好处是代码简单。
根据傅里叶变换的定义式:
X
(
j
ω
)
=
∫
−
∞
∞
x
(
t
)
e
−
j
ω
t
d
t
X(j\omega)=\int_{-\infin}^\infin x(t)e^{-j\omega t}\mathrm{d}t
X(jω)=∫−∞∞x(t)e−jωtdt
可知,连续信号上每个频率傅里叶变换的值就是对信号与指数的乘积做积分的结果。由此,我们可以直接调用matlab的积分函数integral
求值。
首先定义被积函数x=@(t) sin(10.*pi.*t)./(pi.*t).*sin(10.*pi.*t)./(pi.*t).*exp(-1i.*omega.*t);
,之后直接调用积分函数,给定上下限即可。需要注意,虽然matlab的积分函数支持从
−
∞
-\infin
−∞积到
+
∞
+\infin
+∞,但对于许多函数来说,函数值在自变量较大时已经非常小了,对积分的影响微乎其微。对于这类函数我们可以适当缩小积分区间,使得程序运行速度更快。
function res=calc(omega)
x=@(t) sin(10.*pi.*t)./(pi.*t).*sin(10.*pi.*t)./(pi.*t).*exp(-1i.*omega.*t);
res=integral(x,-100,100);
有了上面的函数,我们就能计算出单点的傅里叶变换值,那么只要多求几个点的傅里叶变换,再用平滑曲线连接,我们就可以得到一个较好的频域图像了。
clc;
clear;
n=1001;
lb=-60*pi;
rb=60*pi;
omega=lb:(rb-lb)/(n-1):rb;
X=zeros(1,n);
for i=1:n
X(i)=calc(omega(i));
end
plot(omega,abs(X));
对sinc函数的平方求傅里叶变换得到三角形频域:
对sinc函数求傅里叶变换得到矩形频域: