Issues Translating Custom Discrete Fourier Transform from MATLAB to Python

Question

I'm developing Python software for someone and they specifically requested that I use their DFT function, written in MATLAB, in my program. My translation is just plain not working, tested with sin(2*pi*r). The MATLAB function below:

function X=dft(t,x,f)
% Compute DFT (Discrete Fourier Transform) at frequencies given
%   in f, given samples x taken at times t:
%     X(f) = sum { x(k) * e**(2*pi*j*t(k)*f) }
%             k

shape = size(f);
t = t(:); % Format 't' into a column vector
x = x(:); % Format 'x' into a column vector
f = f(:); % Format 'f' into a column vector

W = exp(-2*pi*j * f*t');
X = W * x;
X = reshape(X,shape);

And my Python interpretation:

def dft(t, x, f):
    i = 1j  #might not have to set it to a variable but better safe than sorry!
    w1 = f * t
    w2 = -2 * math.pi * i
    W = exp(w1 * w2)
    newArr = W * x
    return newArr

Why am I having issues? The MATLAB code works fine but the Python translation outputs a weird increasing sine curve instead of a Fourier transform. I get the feeling Python is handling the calculations slightly differently but I don't know why or how to fix this.

In yout matlab code, you have a `f*t'` which I assume means t transposed. Calculate your `w1` as the sum of what you have at the moment, and see whether that works. -- No, wait. That can't be it. — chw21, Apr 13 '15 at 14:11
@chw21 is right, this is a kronecker produkt between `f` and the transposed `t`. Does Python possibly do an element-wise multiplication there? — hbaderts, Apr 13 '15 at 14:19
Yes, multiplying two numpy arrays returns an array with each element being the product of the two corresponding elements. — chw21, Apr 13 '15 at 14:21
On the difference between numpy arrays and matrices: http://stackoverflow.com/questions/4151128/what-are-the-differences-between-numpy-arrays-and-matrices-which-one-should-i-u/15406440#15406440 — Lee, Apr 13 '15 at 14:27
Works fine for me, you just need to take care of matrix multiplications with `np.dot` as explained by @atomh33ls and also the fact that python uses row major indexing instead of column major indexing in MATLAB. You should be fine. — Divakar, Apr 13 '15 at 14:46
@Divakar the function works for you? I've replaced all the multiplications with np.dot() and did the t.T, but all I'm getting now is the original sine curve (when I plot it, at least). — Trilbador, Apr 13 '15 at 15:03
Did you do the conversion of row to column major indexing at the start of the numpy code with `f = f.T`, `x = x.T` & `t = t.T`? — Divakar, Apr 13 '15 at 15:05
Yes. Though strangely there seems to be no difference between putting that at the top and not having it at all. — Trilbador, Apr 13 '15 at 15:09
Could you list some sample input data: f,x,t so that we could double-check? — Divakar, Apr 13 '15 at 15:56
Is `f` 1d or 2? In MATLAB everything is 2d or higher, in `numpy` arrays can be 1d. Transpose of a 1d does nothing. — hpaulj, Apr 13 '15 at 16:09
t = np.arange(0, 10, 0.005) x = np.sin(2*np.pi*t) f = createArrOfSize(0.0225, 30, size(t)) is what I've been using, where createArrOfSize just creates an array of variable size with the same increment starting with min and ending with max. I'm just getting the original sine curve back from my DFT. — Trilbador, Apr 13 '15 at 23:41

score 1 · Answer 1 · edited May 23 '17 at 12:14

Numpy arrays do element wise multiplication with *.

You need np.dot(w1,w2) for matrix multiplication using numpy arrays (not the case for numpy matrices)

Make sure you are clear on the distinction between Numpy arrays and matrices. There is a good help page "Numpy for Matlab Users":

http://wiki.scipy.org/NumPy_for_Matlab_Users

Doesn't appear to be working at present so here is a temporary link.

Also, use t.T to transpose a numpy array called t.

score 1 · Accepted Answer · answered Apr 14 '15 at 04:56

Here's your MATLAB code -

t = 0:0.005:10-0.005;
x = sin(2*pi*t);
f = 30*(rand(size(t))+0.225);

shape = size(f);
t = t(:); % Format 't' into a column vector
x = x(:); % Format 'x' into a column vector
f = f(:); % Format 'f' into a column vector

W = exp(-2*pi*1j * f*t');  %//'
X = W * x;
X = reshape(X,shape);

figure,plot(f,X,'ro')

And here's one version of numpy ported code might look like -

import numpy as np
from numpy import math
import matplotlib.pyplot as plt

t = np.arange(0, 10, 0.005) 
x = np.sin(2*np.pi*t) 
f = 30*(np.random.rand(t.size)+0.225)

N = t.size

i = 1j
W = np.exp((-2 * math.pi * i)*np.dot(f.reshape(N,1),t.reshape(1,N)))
X = np.dot(W,x.reshape(N,1))
out = X.reshape(f.shape).T

plt.plot(f, out, 'ro')

MATLAB Plot -

enter image description here

Numpy Plot -

enter image description here

Issues Translating Custom Discrete Fourier Transform from MATLAB to Python

2 Answers2