从Matlab的帮助文档中，可以提取出这样一个可以用于实际的快速fourier变换程序。

function [Yf,f] = Spectrum_Calc(yt,Fs)
%功能：实现快速fourier变换
%输入参数：yt为时域信号序列，Fs为采样频率
%返回值：Yf为经过fft的频域序列，f为相应频率

L = length(yt);
NFFT = 2^nextpow2(L);
Yf = fft(yt,NFFT)/L;

Yf = 2*abs(Yf(1:NFFT/2+1));
f = Fs/2 * linspace(0,1,NFFT/2+1);
end

调用方法：

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t));     % Sinusoids plus noise

[Yf,f] = Spectrum_Calc(y,Fs);
plot(f,Yf);

http://s11/middle/84024a4a4cb1499b2200a&690

备注：

（1）关于Fs/2的说明：这个其实是采样定理的规定，即当采样频率fs.max大于信号中最高频率fmax的2倍时(fs.max>=2fmax)，采样之后的数字信号完整地保留了原始信号中的信息，所以，这里f的最大值为Fs/2。如果大于此值，比如取Fs，那么频域上的图像会关于Fs/2对称，而此时靠右侧的尖峰频率实际上是不存在的。

（2）这个Spectrum_Calc是用矩形窗截断的数据，在某些情况下可以加Hanning窗或其他窗，方法很简单

function [Yf,f] = Spectrum_Calc(yt,Fs)
%功能：实现快速fourier变换
%输入参数：yt为时域信号序列，Fs为采样频率
%返回值：Yf为经过fft的频域序列，f为相应频率

L = length(yt);

% 加Hanning窗

w = hann(L);

yt = yt(:).*w(:);

NFFT = 2^nextpow2(L);
Yf = fft(yt,NFFT)/L;

Yf = 2*abs(Yf(1:NFFT/2+1));
f = Fs/2 * linspace(0,1,NFFT/2+1);
end

Matlab帮助文档关于fft的说明

fft

Fast Fourier transform
Syntax

Y = fft(x)
Y = fft(X,n)
Y = fft(X,[],dim)
Y = fft(X,n,dim)

Definitions

The functions Y = fft(x) and y = ifft(X) implement the transform and inverse transform pair given for vectors of length N by:

where

is an Nth root of unity.
Description

Y = fft(x) returns the discrete Fourier transform (DFT) of vector x, computed with a fast Fourier transform (FFT) algorithm.

If the input X is a matrix, Y = fft(X) returns the Fourier transform of each column of the matrix.

If the input X is a multidimensional array, fft operates on the first nonsingleton dimension.

Y = fft(X,n) returns the n-point DFT. fft(X) is equivalent to fft(X, n) where n is the size of X in the first nonsingleton dimension. If the length of X is less than n, X is padded with trailing zeros to length n. If the length of X is greater than n, the sequence X is truncated. When X is a matrix, the length of the columns are adjusted in the same manner.

Y = fft(X,[],dim) and Y = fft(X,n,dim) applies the FFT operation across the dimension dim.

Examples

A common use of Fourier transforms is to find the frequency components of a signal buried in a noisy time domain signal. Consider data sampled at 1000 Hz. Form a signal containing a 50 Hz sinusoid of amplitude 0.7 and 120 Hz sinusoid of amplitude 1 and corrupt it with some zero-mean random noise:
Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t));     % Sinusoids plus noise
plot(Fs*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')

http://s15/middle/84024a4a4cb1499b256ae&690

It is difficult to identify the frequency components by looking at the original signal. Converting to the frequency domain, the discrete Fourier transform of the noisy signal y is found by taking the fast Fourier transform (FFT):
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

http://s5/middle/84024a4a4cb1499c267c4&690

The main reason the amplitudes are not exactly at 0.7 and 1 is because of the noise. Several executions of this code (including recomputation of y) will produce different approximations to 0.7 and 1. The other reason is that you have a finite length signal. Increasing L from 1000 to 10000 in the example above will produce much better approximations on average.