How to produce frequency vector corresponding to FFT() samples?

240 views (last 30 days)
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
  2 Comments
Samuel O'Neill
Samuel O'Neill on 17 May 2021
Update: I'm not sure if my question is clear enough. To be clear: I'm looking for a way to implement the second code block into the first one, so that the function can also produce frequencies of DFT samples. But I'm not sure how to deal with 't' so that it can be implemented into any input
dpb
dpb on 17 May 2021
Edited: 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.

Sign in to comment.

Accepted Answer

David Goodmanson
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.
  2 Comments
dpb
dpb on 17 May 2021
Even though you want to do this on your own w/o FFT(), I suggest reading the example carefully at doc fft that illustrates doing the one-sided PSD via FFT with and without a noise signal and also illustrates one form of normalization to match the time series amplitude.

Sign in to comment.

More Answers (1)

Matt J
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
  2 Comments
Samuel O'Neill
Samuel O'Neill on 17 May 2021
I'm not sure what this does, could you add some comments to explain what you're doing?
Matt J
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.

Sign in to comment.

Categories

Find more on Fourier Analysis and Filtering in Help Center and File Exchange

Products

Community Treasure Hunt

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

Start Hunting!