# Doing DFT without using FFT function

236 views (last 30 days)

Show older comments

Hi all,

I'm working on a project that handles ECG data from arduino and ran into some problems while computing the discrete fourier transform of the ECG. I would like to view the transforms and data collection in real time. While the real time data collection works fine, I would prefer not to use the fft function because for academic uses, the hard coded formula of the fourier transform has more learning value. The code in entirety is as shown below:

clear

delete(instrfindall);

%initialize communication with arduino

arduino = serial('COM6');

arduino.Baudrate=9600;

fopen(arduino);

%initialize figures

figure;

%fig1 = figure;

%fig2 = figure;

%fig3 = figure;

%fig4 = figure;

%initialize sample and counter

sampleSize = 500;

a=1;

while (a < sampleSize) %doubles as a counter, 1000 samples takes about 10seconds

%collects raw data in string format

timeRaw = fgets(arduino);

ecgRaw = fgets(arduino);

%converts string data to numbers

time = str2double(timeRaw) / 1000;

ecg = str2double(ecgRaw);

%data allocation

dataTime(a,1) = time;

dataECG(a,1) = ecg;

%FFT

frequency = 1/time;

fftEcg = fft(ecg);

for k = 1:a

X(k,1) = 0;

for n = 1:a

X(k,1) = X(k,1)+(dataECG(n,1)*exp((-1j)*2*pi*(n-1)*(k-1)/a));

end

end

mag(a,1) = abs(X(a,1));

dataFreq(a,1) = frequency;

dataFft(a,1) = fftEcg;

%low pass filter, set at 50Hz

%sampling freq 80~90Hz, partly determined by delay in arduino

if (frequency < 10)

Y = fftEcg;

dataY(a,1) = Y;

else

Y = 0;

dataY(a,1) = Y;

end

%inverse FFT

dataIfft(a,1) = ifft(Y);

%plotting of figures, subplots are too small, 4 windows too messy

subplot(411);

%figure(fig1);

plot(dataTime,dataECG,'k-');

xlabel('Time (sec)');

ylabel('x(t) magnitude');

subplot(412);

plot(dataFreq,mag,'k-');

xlabel('Frequency (Hz)');

ylabel('X(jw) magnitude');

%subplot(412);

subplot(413);

%figure(fig2);

plot(dataFreq,dataFft,'k-');

xlabel('Frequency (Hz)');

ylabel('X(jw) magnitude');

%subplot(413);

%figure(fig3);

%plot(dataFreq,dataY,'k-');

%xlabel('Frequency (Hz)');

%ylabel('Filtered H(jw) magnitude');

subplot(414);

%figure(fig4);

plot(dataTime,dataIfft,'k-');

xlabel('Frequency (Hz)');

ylabel('Filtered x(t) magnitude');

%command to draw and take in next sample by increasing counter

drawnow;

a = a + 1;

end

%data collection into csv format

%intialize array for final data set

%finalData = zeros(sampleSize,2);

%data collection for raw data only

%for i=1:(sampleSize-1)

% finalData(a,1) = dataTime(a,1);

% finalData(a,2) = dataECG(a,1);

%end

%csvwrite('data.csv', finalData);

fclose(arduino);

return;

In particular the formula that I keyed in is found in this few lines of code:

for k = 1:a

X(k,1) = 0;

for n = 1:a

X(k,1) = X(k,1)+(dataECG(n,1)*exp((-1j)*2*pi*(n-1)*(k-1)/a));

end

end

mag(a,1) = abs(X(a,1));

The results between the fft function and the formula I've input are very different, any ideas how I can modify the formula?

##### 2 Comments

Angus Keatinge
on 28 Apr 2018

Edited: Angus Keatinge
on 28 Apr 2018

This is wrong, the dft is from 0 to N-1 whereas linspace includes the extremities. You won't only have a redundant value at the last index, but every frequency term will be scaled differently to the N point dft (they will be scaled to an N-1 point dft). It is quite a common error, you can correct it by changing the lines:

w = 2*pi*linspace(0,1-1/len,len);

Xk = exp(-1j*w'*(n-1))*xn';

### Accepted Answer

David Young
on 15 Dec 2014

Edited: David Young
on 15 Dec 2014

First, let's confirm that the code you have used for the DFT is correct. Simplifying it a little for clarity (the second subscripts are unnecessary for vectors), we can try it on some test data like this:

N = 20; % length of test data vector

data = rand(N, 1); % test data

X = zeros(N,1); % pre-allocate result

for k = 1:N

X(k) = 0;

for n = 1:N

X(k) = X(k)+(data(n)*exp((-1j)*2*pi*(n-1)*(k-1)/N));

end

end

max(abs(X - fft(data))) % how different from built-in FFT?

This will print out a value of order 10^-14 - i.e. the processes are effectively the same, and this simple DFT code works fine.

So what is wrong is something to do with the way you are using it in within your program. I can see a couple of things that are probably incorrect. One is that a is being used as a loop counter and also as the loop limit in the DFT code. The loop limit needs to be the length of the data vector (which is why I've used N instead of a in my example). The second thing is that you assign the results to a variable called X, but you don't seem to make use of that variable anywhere later on - the results are in effect ignored. fftECG is used, but that isn't given a value anywhere.

So inconsistent use of variables looks like a problem here - remember that each variable that you use needs to be given a value. I think it's very difficult to get code right if you cut and paste snippets from elsewhere.

Finally, note that using the simple DFT code will be far far slower than calling fft() for any significant amount of data.

##### 7 Comments

David Young
on 17 Dec 2014

### More Answers (1)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!