39 views (last 30 days)

Show older comments

I have this function which I wrote to implement the DFT formula without MATLAB's built in fft() function:

function [X] = DFT(x)

N = length(x);

for k = 0:(N-1)

sigma = 0;

for n = 0:(N-1)

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

sigma = sigma + currentVal;

end

X(k+1) = sigma;

end

end

This produces a vector (X) of calculated DFT samples. I now need to modify this function so that I can also produce a vector of the corresponding frequencies. So the function should take x (the input signal) and fs(sampling frequency) as inputs, and return two vectors of the same length: the DFT samples, and their corresponding frequencies. I'm not really sure how I can do this.

As far as I can tell, the code below (which I found on this forum in answer to another post) does what I want. The problem is that I can't see how to implement this into my function so that it can be applied to any input signal. Because in this code, I don't understand what 't' is doing? I've added comments to the code to help explain what I mean.

What is the purpose of 't' here? And how can I modify my function to return frequencies like this code does, when given any input vector 'x'?

Fs = 1000; % sampling rate of 1000 Hz

t = 0:1/Fs:1-1/Fs; %creates a 1x1000 vector from 0 to (1-1/Fs)

x = cos(2*pi*100*t)+randn(size(t)); %confused about this line? what does 't' do? how can I implement this with any input 'x'?

xdft = fft(x); %compute dft of 'x'

xdft = xdft(1:length(x)/2+1); %this line seems to half the vector 'xdft'

DF = Fs/length(x); %this is the increment for the frequency vector

freqvec = 0:DF:Fs/2; %create the frequecy vector

plot(freqvec,abs(xdft)) %plot frequecies against dft output

dpb
on 17 May 2021

function [X] = DFT(x,Fs)

N = length(x);

f = Fs*(0:(N/2))/N;

...

If you don't preallocate X, you'll find it gets really slo0w as N goes up...

Also, length() is dangerous although will work here more often than not, probably. length(x) is defined as max(size(x)) so will return whichever dimension of a 2D array that is the longer, not necessarily the number of rows. If you restrict your usage to vectors only, you're ok; just a caution in general.

David Goodmanson
on 17 May 2021

Edited: David Goodmanson
on 17 May 2021

Hi Samuel,

In the code you borrowed, t is just a time array so that you can create a function of time. The t array starts at 0 (important for the fft) and has array spacing delta_t, By the meaning of Fs as the sampling frequency, delta_t = 1/Fs. Suppose there is a time array with N points and a corresponding time domain function x(t) with N points. Then

t = (0:N-1)*delta_t

= (0:N-1)/Fs

The frequency array also starts at 0 and has array spacing delta_f. You can find delta_f from the golden rule for the fft, which is

delta_t*delta_f = 1/N

So the frequency array is

f = (0:N-1)*delta_f

= (0:N-1)/(N*delta_t)

= (0:N-1)*Fs/N

Usually Fs is known, but if not you can find delta_t directly from the time array. If all you are given is the x array with no accompnying time array and no specified Fs or delta_t, then you can't find the frequency array.

The frequency array and fft output here is full size. One consequence is that you can use ifft to get back to the time domain.

For display purposes you can toss out information by cutting the frequency and output arrays down to (1:N/2+1) as is done in the borrowed code. What remains is zero frequency, all positive frequencies and also the nyquist frequency as the last element. Usually you would take the abs of the fft output and double all the values except for 0 frequency and nyquist frequency. That was not done. The whole cut-in-half procedure is rife with opportunities for error.

dpb
on 17 May 2021

Matt J
on 17 May 2021

Edited: Matt J
on 17 May 2021

[X,freqvec] = DFT(x,Fs);

plot(freqvec,abs(X))

function [X,freqvec] = DFT(x,Fs)

N = length(x);

for k = 0:(N-1)

sigma = 0;

for n = 0:(N-1)

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

sigma = sigma + currentVal;

end

X(k+1) = sigma;

end

X=iffshift(X);

if nargout>1

freqvec=( (0:N-1)-ceil((N-1)/2) )/N*Fs;

end

end

Matt J
on 17 May 2021

I've just added a few lines to your posted code, mainly the part that produces the frequency vector.

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

Start Hunting!