Matlab 版 (精华区)

发信人: candle ( 马 走 日), 信区: Matlab
标  题: How can I use FFT to obtain simple spectral analy
发信站: 紫 丁 香 (Fri Dec 24 17:48:36 1999), 转信

1702
              How Can I Use FFT to Obtain Simple
              Spectral Analysis Plots?
    Revison: 1.1
                                           Last Date Modified: 29-June-1999
    Question:
    How can you correctly scale the output of the FFT function to
    obtain a meaningful magnitude versus frequency plot?
    Answer
    Assume "x" is a vector containing your data. A sample vector used
    with this technical note is
        Fs=1000;                  % sampling frequency
        Fn=Fs/2;                  % Nyquist frequency
        t=0:1/Fs:1;               % time vector sampled at Fs Hz,
                                  % length of 1 second
        x = sin(2*pi*t*200);      % sine wave of 200 Hz.
    First, you need to call the fft function. For the fastest possible
    ffts, you will want to pad your data with enough zeros to make its
    length a power of two. The fft function does this for you
    automatically if you give a second argument specifying the overall
    length of the fft, as demonstrated below:
        % Next highest power of 2 greater than or equal to
        % length(x):
        NFFT=2.^(ceil(log(length(x))/log(2)));
        % Note: If you are using version 3.0 or above of the
        % Signal Processing Toolbox, use:
        %
        %         NFFT = 2^(nextpow2(length(x)))
        % Take fft, padding with zeros, length(FFTX)==NFFT
        FFTX=fft(x,NFFT);.
    If NFFT is even (which it will be, if the above two commands are
    used), then magnitude of the fft will be symmetric such that the
    first (1+NFFT/2) points are unique, and the rest are symmetrically
    redundant. FFTX(1) is the DC component of x, and
    fftx(1+NFFT/2) is the Nyquist frequency component of x.If NFFT
    is odd, however, the Nyquist frequency component is not
    evaluated, and the number of unique points is (NFFT+1)/2. This
    can be generalized for both cases to ceil((NFFT+1)/2).
        NumUniquePts = ceil((NFFT+1)/2);
        % fft is symmetric, throw away second half
        FFTX=FFTX(1:NumUniquePts);
    Next, we calculate the magnitude of the fft:
        MX=abs(FFTX);            % Take magnitude of X
        % Multiply by 2 to take into account the fact that we
        % threw out second half of FFTX above
        MX=MX*2;
    Note that simply multiplying MX by 2 here is the quick-and-dirty
    solution. Technically, the DC component and the Nyquist
    component (if it exists, see above) are unique and should not be
    multiplied by 2. To be accurate, the following lines can be added:
        MX(1)=MX(1)/2;
        if ~rem(NFFT,2),
             MX(length(MX))=MX(length(MX))/2;
        end
    Next, we need to take into account the fact that MATLAB doesn't
    scale the output of fft by the length of the input:
        % Scale the FFT so that it is not a function of the
        % length of x
        MX=MX/length(x);
    Now we create a frequency vector:
        % This is an evenly spaced frequency vector with
        % NumUniquePts points.
        f=(0:NumUniquePts-1)*2/NFFT;
        % Multiply this by the Nyquist frequency
        % (Fn ==1/2 sample freq.)
        f=f*Fn;
    The above two lines can be combined to
        f=(0:NumUniquePts-1)*2*Fn/NFFT;
    Finally, we can generate the plot
        plot(f,MX);
    Bringing this all together, we get the following M-file:
        Fs=1000;                  % sampling frequency
        Fn=Fs/2;                  % Nyquist frequency
        t=0:1/Fs:1;               % time vector sampled at Fs Hz,
                                  % length of 1 second
        x = sin(2*pi*t*200);      % sine wave of 200 Hz.
        % Next highest power of 2 greater than or equal to
        % length(x):
        NFFT=2.^(ceil(log(length(x))/log(2)));
        % Take fft, padding with zeros, length(FFTX)==NFFT
        FFTX=fft(x,NFFT);
        NumUniquePts = ceil((NFFT+1)/2);
        % fft is symmetric, throw away second half
        FFTX=FFTX(1:NumUniquePts);
        MX=abs(FFTX);            % Take magnitude of X
        % Multiply by 2 to take into account the fact that we
        % threw out second half of FFTX above
        MX=MX*2;
        MX(1)=MX(1)/2;   % Account for endpoint uniqueness
        MX(length(MX))=MX(length(MX))/2;  % We know NFFT is even
        % Scale the FFT so that it is not a function of the
        % length of x.
        MX=MX/length(x);                  %
        f=(0:NumUniquePts-1)*2*Fn/NFFT;
        plot(f,MX);


--
---------------------------------------------------
I BELIEVE I CAN FLY. I BELIEVE I CAN TOUCH THE SKY.
-------------------   CANDLE   --------------------

※ 来源:．紫 丁 香 bbs.hit.edu.cn．[FROM: 150.59.34.186]