不难，有空可以给你看看

真得很感谢，希望能和您成为朋友，已经加了QQ号了
是用MATLAB编，导师告诉我A=(a(f).^2+b(f).^2).^1/2,我编完后怎么也和书上不一样
程序如下：1。f=603频谱
clear all
clc
Af1=[];Bf1=[];
for f=601:4/82:605;
    k1=0;
    for i=0:1023;
        x1=double(1000*cos(2*pi*124.5*i)+2000*cos(2*pi*603*i+0.4*pi));
        k1=k1+x1*cos(2*pi*i*f/2048);
    end
    Af1=[Af1,2/1024*k1];
     f
     k2=0;
    for i=0:1023;
        x2=double(1000*cos(2*pi*124.5*i)+2000*cos(2*pi*603*i+0.4*pi));
        k2=k2+x2*sin(2*pi*i*f/2048);
    end
    Bf1=[Bf1,2/1024*k2];
     f   
 end
 A=((Af1).^2+(Bf1).^2).^(1/2)
 figure(5)
 plot([601:4/82:605],abs(A));




2。f=124.5频谱
clear all
clc
Af1=[];Bf1=[];
for f=122:4/82:126;
    k1=0;Af1=[];
    for i=0:1023;
        x1=double(1000*cos(2*pi*124.5*i)+2000*cos(2*pi*603*i+0.4*pi));
        k1=k1+x1*cos(2*pi*i*f/2048);
    end
    Af1=[Af1,2/1024*k1];
     f
     k2=0;
    for i=0:1023;
        x2=double(1000*cos(2*pi*124.5*i)+2000*cos(2*pi*603*i+0.4*pi));
        k2=k2+x2*sin(2*pi*i*f/2048);
    end
    Bf1=[Bf1,2/1024*k2];
     f   
 end
 A=((Af1).^2+(Bf1).^2).^(1/2)
 figure(5)
 plot([122:4/82:126],abs(A));

基本上和书本上一样的。其余的，你参考下

% Author : Chianel.Xie
% Date   : 15/11/2007
function FFT_FS
% Para. Initialize
fs = 2560; % HZ. The Sample Frequency
N  = 1024; % The Sample Number
fc = fs/2.56; % The Stop Frequency
t  = 1/fs:1/fs:N/fs; % The Signal Time Duration
x  = 5*sin(2*pi*64*t) + 5*cos(2*pi*53*t) + 10*sin(2*pi*89*t);% Signal

%%%%%%%%%%--------%%%%%%%%%%%%%%%
%            FFT
%%%%%%%%%%--------%%%%%%%%%%%%%%%
% 1024-FFT
Y = fft(x,N);
% The power spectrum
Pyy = Y.* conj(Y) / N; % OR: Pyy = abs(Y)/N;
% Get part frequency
f = fs*(0:N/2)/N; % OR:

Flag = 1;
% Plot Figure
if Flag == 1
    figure
    plot(f,Pyy(1:N/2+1))
    axis([0 150 0 2e4])
    title('Frequency content of x')
    xlabel('frequency (Hz)')
    grid
end