信息处理 - DFT 的维基百科方程的实现 - 吾爱随笔录

DFT 的维基百科方程的实现

信息处理傅里叶变换 Python 自由度

2021-12-25 11:23:56

我正在写一个简单的傅里叶变换实现，并查看了维基百科上的 DFT 方程以供参考，当我注意到我在做一些不同的事情时，经过思考后觉得维基百科的版本一定是错误的，因为它很容易想到一个信号表明，当傅立叶变换（使用该等式）将返回不正确的频谱：因为该等式仅将信号围绕复平面包裹一次（由于的），任何周期性的信号偶数次（在包裹复平面时）将没有频谱，因为在 DFT 期间出现的通常峰值（在围绕单位圆时）将相互抵消（当它们出现偶数时）。 $n/N$ $0<n<N-1$

为了检查这一点，我编写了一些生成以下图像的代码，这似乎证实了我的想法。

“时间使用方程”使用方程与向量时间（例如采样的它可以在下面的函数中找到。

X_{f} = \sum_{n = 0}^{N - 1} x_{n} (\cos (2 π f t_{n}) - i \sin (2 π f t_{n}))

$X_f = \sum_{n=0}^{N-1} x_n (\cos(2\pi f t_n) - i \sin(2\pi f t_n))$

t

$t$

t_{n}

$t_n$

x_{n}

$x_n$ ft

上面链接的维基百科方程在这里复制以供参考：可以在函数中找到。

X_{f} = \sum_{n = 0}^{N - 1} x_{n} (\cos (2 π f \frac{n}{N}) - i \sin (2 π f \frac{n}{N}))

$X_f = \sum_{n=0}^{N-1} x_n \left(\cos\left(2\pi f \frac{n}{N}\right) - i\sin\left(2\pi f \frac{n}{N}\right)\right)$ ft2

import numpy as np
import matplotlib.pyplot as plt 
plt.style.use('ggplot')

def ft(t, s, fs):
    freq_step = fs / len(s)
    freqs = np.arange(0, fs/2 + freq_step, freq_step)
    S = []
    for freq in freqs:
        real = np.sum(s * np.cos(2*np.pi*freq * t)) 
        compl = np.sum(- s * np.sin(2*np.pi*freq * t)) 
        tmpsum = (real**2 + compl**2) ** 0.5 
        S.append(tmpsum)
    return S, freqs

def ft2(s, fs):  # Using wikipedia equation
    nump=len(s)
    freq_step = fs / nump
    freqs = np.arange(0, fs/2 + freq_step, freq_step)
    S = []
    for i, freq in enumerate(freqs):
        real = np.sum(s * np.cos(2*np.pi*freq * i/nump))
        compl = np.sum(- s * np.sin(2*np.pi*freq * i/nump))
        tmpsum = (real**2 + compl**2) ** 0.5 
        S.append(tmpsum)
    return S, freqs


def main():
    f = 5 
    fs = 100 
    t = np.linspace(0, 2, 200)
    y = np.sin(2*np.pi*f*t) + np.cos(2*np.pi*f*2*t)

    fig = plt.figure()
    ax = fig.add_subplot(311)
    ax.set_title('Signal in time domain')
    ax.set_xlabel('t')
    ax.plot(t, y)

    S, freqs = ft(t, y, fs) 

    ax = fig.add_subplot(312)
    ax.set_xticks(np.arange(0, freqs[-1], 2)) 
    ax.set_title('Time using equation')
    ax.set_xlabel('frequency')
    ax.plot(freqs, S)

    S, freqs = ft2(y, fs) 
    ax = fig.add_subplot(313)
    ax.set_title('Using Wiki equation')
    ax.set_xlabel('frequency')
    ax.set_xticks(np.arange(0, freqs[-1], 2)) 
    ax.plot(freqs, S)

    plt.tight_layout()
    plt.show()

main()

显然，我似乎不太可能在如此高调的 wiki 页面上随机发现错误。但我看不出我所做的有什么错误？

3个回答

你有一个错误ft2。你正在增加i，并freq在一起。这不是你希望你的求和工作的方式。我搞砸了修复它，但它变得一团糟。我决定从离散的角度重写它，而不是尝试使用连续的术语。在 DFT 中，采样率是无关紧要的。重要的是使用了多少样本 ( N)。然后，bin 编号 ( k) 对应于以每帧周期为单位的频率。我试图让您的代码尽可能完整，以便您可以轻松理解它。我还展开了 DFT 计算循环，希望能更好地揭示它们的性质。

希望这可以帮助。

赛德

将 numpy 导入为 np
将 matplotlib.pyplot 导入为 plt

定义 ft(t, s, fs):
    freq_step = fs / len(s)
    频率 = np.arange(0, fs/2, freq_step)
    S = []
    对于频率中的频率：
        实数 = np.sum(s * np.cos(2*np.pi*freq * t))
        compl = np.sum(- s * np.sin(2*np.pi*freq * t))
        tmpsum = (real**2 + compl**2) ** 0.5
        S.append(tmpsum)
    返回 S，频率

def ft3(s, N): # 更有效的维基百科方程形式

    S = []

    切片 = 0.0
    条子 = 2*np.pi/float(N)

    对于范围内的k（N / 2）：

        sum_real = 0.0    
        sum_imag = 0.0
        角度 = 0.0
        对于范围内的n（N）：
            sum_real += s[n] * np.cos(角度)
            sum_imag += -s[n] * np.sin(角度)
            角度 += 切片

        切片 += 条子
        tmpsum = (sum_real**2 + sum_imag**2) ** 0.5
        S.append(tmpsum)

    返回 S

def ft4(s, N): # 使用维基百科方程

    S = []

    对于范围内的k（N / 2）：

        sum_real = 0.0    
        sum_imag = 0.0
        对于范围内的n（N）：
            sum_real += s[n] * np.cos(2*np.pi*k*n/float(N))
            sum_imag += -s[n] * np.sin(2*np.pi*k*n/float(N))

        tmpsum = (sum_real**2 + sum_imag**2) ** 0.5
        S.append(tmpsum)

    返回 S

def ft5(s, N): # 统一加权和的根

    条子 = 2 * np.pi / float( N )

    root_real = np.zeros(N)
    root_imag = np.zeros( N )

    角度 = 0.0
    对于范围内的 r (N)：
        root_real[r] = np.cos( 角度 )
        root_imag[r] = -np.sin( 角度 )
        角度 += 条子

    S = []

    对于范围内的k（N/2）：

        sum_real = 0.0    
        sum_imag = 0.0
        r = 0

        对于范围内的n（N）：
            sum_real += s[n] * root_real[r]
            sum_imag += s[n] * root_imag[r]
            r += k
            如果 r >= N : r -= N

        tmpsum = np.sqrt( sum_real*sum_real + sum_imag*sum_imag )
        S.append(tmpsum)

    返回 S

定义主（）：

    N = 200
    fs = 100.0

    time_step = 1.0 / fs
    t = np.arange(0, N * time_step, time_step)

    f = 5.0
    y = np.sin(2*np.pi*f*t) + np.cos(2*np.pi*f*2*t)

    无花果 = plt.figure()
    斧头 = fig.add_subplot(311)
    ax.set_title('时域信号')
    ax.set_xlabel('t')
    ax.plot(t, y)

    S, 频率 = ft(t, y, fs)

    斧头 = fig.add_subplot(312)
    ax.set_xticks(np.arange(0, freqs[-1], 2))
    ax.set_title('时间使用方程')
    ax.set_xlabel('频率')
    ax.plot（频率，S）

    S = ft3(y, N)
    斧头 = fig.add_subplot(313)
    ax.set_title('使用维基方程')
    ax.set_xlabel('频率')
    ax.set_xticks(np.arange(0, freqs[-1], 2))
    打印 len(S), len(freqs)
    ax.plot（频率，S）

    plt.tight_layout()
    plt.show()

主要的（）

我不会看你的代码。维基百科页面看起来不错，但它是数学家和电气工程师之间“格式战争”或“符号战争”或“风格战争”的一个很好的例子。其中一些，我认为数学家是对的。EE 不应该采用“ ”作为虚数单位。也就是说，这是 DFT 的更好表达，反之则为： $j$

DFT：

X [k] = \sum_{n = 0}^{N - 1} x [n] e^{- j 2 π n k / N}

$X[k] = \sum\limits_{n=0}^{N-1} x[n] \, e^{-j 2 \pi nk/N}$

iDFT:

x [n] = \frac{1}{N} \sum_{k = 0}^{N - 1} X [k] e^{j 2 π n k / N}

$x[n] = \frac{1}{N} \sum\limits_{k=0}^{N-1} X[k] \, e^{j 2 \pi nk/N}$

因为做 DSP 的电气工程师喜欢使用作为“时间”中的样本序列，而作为“频率”中的离散样本序列。数学家可能更喜欢这个： $x[n]$ $X[k]$

DFT：

X_{k} = \sum_{n = 0}^{N - 1} x_{n} e^{- i 2 π n k / N}

$X_k = \sum\limits_{n=0}^{N-1} x_n \, e^{-i 2 \pi nk/N}$

iDFT:

x_{n} = \frac{1}{N} \sum_{k = 0}^{N - 1} X_{k} e^{i 2 π n k / N}

$x_n = \frac{1}{N} \sum\limits_{k=0}^{N-1} X_k \, e^{i 2 \pi nk/N}$

这与维基百科页面相同。

您可能需要更多地注意在指数中使用或以及如何将其转换为或与术语相反。 $+$ $-$ $+$ $-$ $\sin(\cdot)$

我回到这一点并尝试派生有助于使事情更有意义的离散版本：

不知何故 $f_k t_n = f(n, k, N)$

$f_k = \frac{f_s}{N}k$ 和 $t_n = \frac{T}{N}n$

$f_s = \frac{N}{T}$

所以

$f_k t_n = \frac{f_s}{N}k\frac{T}{N}n = \frac{N}{TN}k\frac{T}{N}n = \frac{kn}{N}$

完毕！

其它你可能感兴趣的问题

上一篇过滤流式音频下一篇如何将函数窗口应用于信号？