Doing DFT without using FFT function

Question

Ben 2014-12-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/166660-doing-dft-without-using-fft-function

回答： Piotr Gregor 2022-7-1

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 个评论
显示无隐藏无

Soroush 2016-2-15

编辑：Soroush 2016-2-15

在 MATLAB Online 中打开

One Dimensional DFT Function, without using for loop:

function [Xk] = dft(xn)
len = length(xn);
w = 2*pi*linspace(0,1,len);
n = 1:len;
Xk = exp(-1j*w'*n)*xn';

Angus Keatinge 2018-4-28

编辑：Angus Keatinge 2018-4-28

在 MATLAB Online 中打开

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';

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

David Young 2014-12-15

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/166660-doing-dft-without-using-fft-function#answer_162383

编辑：David Young 2014-12-15

在 MATLAB Online 中打开

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 个评论
显示 5更早的评论隐藏 5更早的评论

Ben 2014-12-15

在 MATLAB Online 中打开

I would agree on what you have said mathematically but I used the built in fft function on one incoming value from the serial port and it managed to compute fft in real time, because the figures are able to be updated live.

This is where I used the fft function after collecting data from the serial port 'arduino':

        %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);

and also

dataFft(a,1) = fftEcg;

Output onto a graph is done here:

          %subplot(412);
          subplot(413);
          %figure(fig2);
          plot(dataFreq,dataFft,'k-');
          xlabel('Frequency (Hz)');
          ylabel('X(jw) magnitude');

But however when I try to represent this built in function with a self-written DFT code as per discussed earlier, the calculated values in the array X turns out totally different from the values calculated by the inbuilt fft() function.

Which is the reason why I questioned if the equation I've input is wrong.. Or is there is something else that the inbuilt fft() does besides calculating DFT?

Ben 2014-12-16

在 MATLAB Online 中打开

Hm, it is strange because ecg is just a double. It isn't even an array so there isn't a length to speak of. Which hence prompted me to believe fft() does not actually apply to a fixed vector.

On a side note, I had to use a in each of the loops because in computation of DFT, a is synonymous to the sample size N, which you had pointed out correctly should be defined beforehand conventionally. However, because the size of dataFft is always increasing as the serial port receives an additional data point, N (and therefore a) will always be increasing as well if I were to carry out real-time data collection instead of declaring the length beforehand.

I am puzzled because the fft() function was able to work this out on its own without me declaring a fixed length for it.

This was the code BEFORE I implemented the DFT hard-coded formula and worked perfectly fine:

clear
delete(instrfindall);
%initialize communication with arduino
arduino = serial('COM6');
arduino.Baudrate=9600;
fopen(arduino);
%initialize figures
figure;
%initialize sample and counter
sampleSize = 1000;
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);
      %FFT
      frequency = 1/time;
      fftEcg = fft(ecg);
      %data allocation
      dataTime(a,1) = time;
      dataECG(a,1) = ecg;
      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);
      plot(dataTime,dataECG,'k-');
      xlabel('Time (sec)');
      ylabel('x(t) magnitude');
      subplot(412);
      plot(dataFreq,dataFft,'k-');
      xlabel('Frequency (Hz)');
      ylabel('X(jw) magnitude');
      subplot(413);
      plot(dataFreq,dataY,'k-');
      xlabel('Frequency (Hz)');
      ylabel('Filtered H(jw) magnitude');
      subplot(414);
      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(i,1) = dataTime(i,1);
%    finalData(i,2) = dataECG(i,1);
%end
%csvwrite('data.csv', finalData);
fclose(arduino);
return;

And because fft() is a shortcut function and has no academic value whatsoever, a hardcoded formula is preferred for teaching purposes, despite the fact that fft() will be computationally faster than any form of hardcoded DFTs.

Ben 2014-12-17

ah, I found out where it went wrong. ecg should be stored as an array and not be taken as a singular value.

Many thanks!

David Young 2014-12-17

That sounds right. fft() will work on a single value (a scalar) - it just treats it as a vector of length 1, and returns its argument. This is correct mathematically, since the DFT of a scalar is just the scalar itself.

请先登录，再进行评论。