在不使用内置函数的情况下在 python 上编写离散傅立叶变换

Question

这是我的东西

import math
import numpy as np
from cmath import sqrt

# Ask the user for the number of elements in the sequence
n = int(input("n: "))

# create a matrix full of n elements
x = [0]*n
for i in range(0,n):
    x[i] = int(input(f"x{i}: "))

j= sqrt(-1)

y = [0]*n
for k in range(0,n):
    for l in range(0,n):
        y[k] += x[l]*(np.exp((-np.pi *j* k*2)/n))

# Print out the answers
for i in range(0,n):
print(f"y{i} = {y[i]}")

例一

我使用

x = {0,1,2,3}

得到的答案是：

y0 = (6+0j)
y1 = (3.6739403974420594e-16-6j) WRONG
y2 = (-6-7.347880794884119e-16j) WRONG
y3 = (-1.102182119232618e-15+6j) WRONG

应该是：

6, -2+2j, -2, -2-2j

在这里部分派生：

示例 2

如果 n 是 2 并且

x = {1,2}

那么预期的答案是：

3/sqrt(2)

和

-1/sqrt(2)

示例 3`

如果 n 是 4 并且

x = {1,1,0,0}

那么答案是

2

、

1-j

、

0

和

1+j

---

这里有一个修复，但是当我做

x={1,2}

部分时，我没有得到理论上的答案......它们非常不同

import math
import numpy as np
from cmath import sqrt
n = int(input("n: "))

x = [0]*n

#Create a matrix for n amount of inputs
for i in range(0,n):
    x[i] = int(input(f"x{i}: "))

#Empty matrix for the answers for later

y = [0]*n
# define j so that it can be used later
j = sqrt(-1)
for k in range(0,n):
    for l in range(0,n):
        # Carry out the summation
        y[k] += x[l]*(np.exp((-np.pi *j*l* k*2)/n))

for i in range(0,n):
    print(f"y{i} = {y[i]}")

Answer 1

更快的方法是使用矩阵乘积（实际上是一个 numpy 数组和

matmul

函数）：

# Custom matrix
import numpy as np
k = np.arange(N)
M = np.exp(-2j * np.pi * k[:, None] * k / N)
X = np.matmul(xn, M)

xn

是信号样本，N是样本数。

作为参考，有一个返回相同数组的 scipy 函数：

# Built-in matrix
import scipy.linalg as sl
M = sl.dft(N)
X = np.matmul(xn, M)

这段代码显示了这两种方法与 FFT 的等价性：

import numpy as np
import matplotlib.pyplot as plt

def show_peak_frequencies(ax, ff):
    lims = ax.get_xlim()
    tax = ax.twiny()
    ticks = np.hstack([-ff, ff])
    tax.set_xticks(ticks)
    tax.set_xlim(lims)
    tax.grid()

mag = lambda X: abs(sf.fftshift(X))
omega = lambda f: 2*np.pi * f

# Sampling parameters
N = 70
fs = 10
Ts = 1 / fs
max_t = (N-1) * Ts

# Signal components
ff = np.array([1, 2.5, 3.5]) # frequencies Hz
w1, w2, w3 = omega(ff)

# Continuous signal
t = np.linspace(0, max_t, 200)
x = lambda t: 2 * np.cos(w1*t) + np.cos(w2*t) + 0.5 * np.cos(w3*t)
xt = x(t)

# Actual samples
n = np.arange(N)
xn = x(n*Ts)

fig, axes = plt.subplots(ncols=2, squeeze=True, figsize=(9, 3), layout='constrained')

ax = axes[0]
ax.plot(t, xt, c='gray')
ax.stem(n*Ts, xn)
ax.set_title(f'Signal sampled at {1/Ts:0.2f}Hz')
ax.set_xlabel('t (s)')

# Built-in fft
import scipy.fft as sf
X1 = sf.fft(xn)
f = sf.fftfreq(N, Ts)
f = sf.fftshift(f)

# Built-in matrix
import scipy.linalg as sl
M2 = sl.dft(N)
X2 = np.matmul(xn, M2)

# Custom matrix
k = np.arange(N)
M3 = np.exp(-2j * np.pi * k[:, None] * k / N)
X3 = np.matmul(xn, M3)

ax = axes[1]
ax.plot(f, mag(X1), label='FFT', c='lime', lw=7, alpha=0.5)
ax.plot(f, mag(X2), label='scipy matrix', c='navy', lw=4, ls='--')
ax.plot(f, mag(X3), label='custom matrix', c='magenta')
show_peak_frequencies(ax, ff)
ax.set_xlabel('f (Hz)')
ax.legend()