MATLAB Answers

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

조회 수: 2(최근 30일)
Samuel O'Neill
Samuel O'Neill 2021년 5월 17일
편집: dpb 2021년 5월 17일
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
dpb
dpb 2021년 5월 17일
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
David Goodmanson 2021년 5월 17일
편집: David Goodmanson 2021년 5월 17일
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
dpb
dpb 2021년 5월 17일
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.

댓글을 달려면 로그인하십시오.

추가 답변(1개)

Matt J
Matt J 2021년 5월 17일
편집: Matt J 2021년 5월 17일
[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
Matt J
Matt J 2021년 5월 17일
I've just added a few lines to your posted code, mainly the part that produces the frequency vector.

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

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

Start Hunting!

Translated by