MATLAB: How to produce frequency vector corresponding to FFT() samples

dftfftfourierfrequencyMATLABsampling

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

Best Answer

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.
Related Question