信息处理 - 用 FFT 逼近洛伦兹傅里叶变换 - 吾爱随笔录

用 FFT 逼近洛伦兹傅里叶变换

信息处理 fft 傅里叶变换

2022-02-04 02:37:25

我有一个洛伦兹频率分布

$F(w) = \frac{1+iz}{1+z^2}$

在哪里

$z = \frac{w-\Omega}{R}$

和 $\Omega$ 是峰值频率，R 是衰减常数。我知道分析傅里叶变换应该是

$F(t) = exp(i\Omega 2\pi t)exp(-Rt)$

当我对该表达式进行 FFT 时，它不会返回原始频率图。我知道可能有一个比例因子（1/n）应该在某个地方，但即使我缩放一个频率，如果我然后改变 $\Omega$ 或 R，幅度不再正确缩放，表明缩放因子是 $\Omega$ 和/或 R。FFT 似乎也沿强度轴镜像。

我对 DSP 很陌生，但我知道连续傅立叶变换不是离散傅立叶变换。我读过这篇文章（https://dspillustrations.com/pages/posts/misc/approximating-the-fourier-transform-with-dft.html），但这种方法使近似值变得更糟。

当我对时间信号进行 FFT 时，我想返回原始频率分布。我错过了一些基本的东西还是一个相当简单的缩放错误？我在下面附上了我的代码。

干杯。

# R script to compare FFT and Analytical fourier transform
library(SynchWave)

#-------------------------------------------------
# Frequency and time axes
n <- 100
f <- seq(0, 1, length.out = n)
t <- seq(0, n, length.out = n)

# peak paramaters 
O <- 0.3 # Frequency values from 0->1
R <- 0.04 # Decay in arbritrary units

z <- (f-O)/R

# The original lorentzian frequency 
ff <-complex(re = 1, im = z)/(1 + z^2)

# creating the time domain signal
ftideal <- exp(-R*t)*exp(complex(i = (O)*2*pi*t))

unscaled <- (fft(ftideal))
scaled <- unscaled - min(Re(unscaled))

plot(f, Re(ff), type = 'l')
lines(f, Re(scaled), type = "l", col = 'red')

```

1个回答

首先，看看这个非常相关的问题及其答案。其次，您的分析解决方案是错误的，因此您没有看到正确的结果也就不足为奇了。

时域函数

\begin{matrix} (1) & x (t) = e^{j 2 π f_{0} t} e^{- R t} \end{matrix}

$x(t)=e^{j2\pi f_0t}e^{-Rt}\tag{1}$

具有以下傅里叶变换：

\begin{matrix} (2) & X (f) = \frac{1}{R} \frac{1 - j \frac{2 π (f - f_{0})}{R}}{1 + {(\frac{2 π (f - f_{0})}{R})}^{2}} \end{matrix}

$X(f)=\frac{1}{R}\frac{1-j\frac{2\pi(f-f_0)}{R}}{1+\left(\frac{2\pi(f-f_0)}{R}\right)^2}\tag{2}$

此外，缩放应该是乘法而不是加法（除非您在对数域中），并且您还应该不仅查看复值函数的实部。

正如这个答案中所解释的，用 DFT 近似 CTFT 通常会导致两种类型的错误：截断误差（由于时域函数的截断）和混叠误差（由于对时域函数进行采样）。通过选择较大的采样频率和较大的 DFT 长度，可以使这些误差变小。

下面是 Octave/Matlab 中的示例代码，显示了给定函数的 CTFT 如何通过 DFT 近似：

F0 = 0.3；
R = 0.04；

Fs = 30; ％ 采样频率
Ts = 1/Fs；
Tmax = 100；% 时域信号长度
N = 轮次（Tmax/Ts）；% DFT 长度

t = (0:N-1) * Tmax / (N-1);
N2 = 圆形（N/2）；
f = (0:N2-1)/N*Fs；

乐趣 = exp(1i*2*pi*F0*t) .* exp(-R*t);

CTFT 的 % 解析表达式
z = 2*pi*(f - F0) / R；
FTfun = (1 - 1i*z) ./ (1 + z.^2) / R;

% DFT 近似 CTFT
FTfun2 = Ts * fft(fun,N);
FTfun2 = FTfun2(1:N2);

情节（f，20 * log10（abs（FTfun）），f，20 * log10（abs（FTfun2）））
    标题（'幅度（dB）'），xlabel（'f'），图例（'CTFT'，'DFT'），网格上
    轴（[0,Fs/2,-40,30]）

其它你可能感兴趣的问题

上一篇使用 FFT 进行串行滤波下一篇使用两个 USRP 发送/接收 BPSK 信号时遇到问题