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

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

[Samuel O'Neill]

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]

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.

Answer 1

[David Goodmanson]

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.

[Samuel O'Neill]

Thankyou thats a super helpful answer :)

[dpb]

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.

Answer 2

[Matt J]

[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

[Samuel O'Neill]

I'm not sure what this does, could you add some comments to explain what you're doing?

[Matt J]

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