Discrete Fourier transforms don’t work/are inefficient in python
I’m trying to write a discrete Fourier transform function in python that takes the spectral density of the signal as an array (which I’ll then output graphically). I’m using matrix multiplication to do this. My code seems to work for a small portion of the data, but it takes a long time to process; For large amounts of data, such as a WAV file, it can never complete the task. The current features are:
from scipy.io import wavfile
import numpy as np
import matplotlib.pyplot as plt
def esd (data, fs):
a=[]
for i in range(int(np.size(data)/100)):
dt = 1/(fs)
fvec = np.arange(100*i , 100+(100*i) , 1)
nvec = np.arange(0 , np.size(data) , 1)
Wconst = np.exp(-1j*2*np.pi/np.size(data))
ematrix = Wconst**(np.outer(fvec,nvec))
b = (dt**2)*(abs(np.inner(ematrix , data))**2)
a = np.concatenate((a,b))
return a
where data is
the data set and fs is the sampling frequency. Is there a reason why this feature is so slow or inefficient? It is designed to process signals in 100 frequency blocks, otherwise the matrix would be very large.
Solution
The algorithm implements the “naïve” discrete Fourier transform (DFT) by calculating the Vandermonde frequency matrix. Here’s the problem:
b = (dt ** 2) * abs(np.inner(ematrix, data)) ** 2
This simple implementation uses direct matrix vector multiplication with a run time of O(N**2), where N == data.size. So, when you get bigger data, it gets worse and may not be done for your WAV file in a reasonable amount of time.
That’s the whole point of using the fast Fourier transform algorithm, which makes use of a large number of structures in the problem. Most notably, FFT recursively breaks down the problem into smaller problems of size N/2. This recursive structure makes the FFT run time O(N log N).