aw_psd.m

This getting started exercise will guide you through the step-by-step process of transforming a MATL

%       AccelDSP 8.1.1 build 690 Production, compiled Apr 26 2006 
% 
%    THIS IS UNPUBLISHED, LICENSED SOFTWARE THAT IS THE CONFIDENTIAL 
%        AND PROPRIETARY PROPERTY OF XILINX OR ITS LICENSORS 
% 
%     Copyright(c) Xilinx, Inc., 2000-2006, All Rights Reserved. 
%     Reproduction or reuse, in any form, without the explicit written 
%     consent of Xilinx, Inc., is strictly prohibited. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Description : Custom PSD function used in FIR top level script.            % 
%%                                                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [P,F]=aw_psd(X,NFFT,Fs)
if (min(size(X))~=1)
    error('aw_psd only operates on row or column vectors')
end
% Convert input to column vector
X=X(:);
% Set default values if not specified
if (nargin < 3)
    Fs = 2;
end
if (nargin < 2)
    if (length(X)<256)
        NFFT = length(X);
    else
        NFFT = 2^8;
    end
end
RealVal = isreal(X);
if (RealVal == 0)
    Fs=Fs*2;
end
% Zero pad input if needed
if (length(X) < NFFT)
    X = [X;zeros(NFFT-length(X),1)];
end
% Determine available averaging
if (length(X) > NFFT)
    AVE = floor(length(X)/NFFT);
    X = X(1:AVE*NFFT);
%    X = [X;zeros(AVE*NFFT-length(X),1)];
else
    AVE = 1;
end
X=reshape(X,NFFT,AVE);
%Create Hanning Window
if ~rem(NFFT,2)
    % Even length window
    m=(NFFT)/2;
    w = .5*(1 - cos(2*pi*(1:m)'/(NFFT+1)));
    w = [w; w(end:-1:1)];
else
    % Odd length window
    m=(NFFT+1)/2;
    w = .5*(1 - cos(2*pi*(1:m)'/(NFFT+1)));
    w = [w; w(end-1:-1:1)];
end
hm=kron(ones(1,AVE),w);
%Compute PSD
P=abs(fft(X.*hm)).^2/((length(X))/2);
P=mean(P,2)*1.332;
if (RealVal == 1)
    P=P(1:end/2,:);
    NFFT=NFFT/2;
end
F=linspace(0,1-1/NFFT,NFFT)*Fs/2;
if (nargout == 0)
    plot(F,10*log10(P),'r')
   grid on
    xlabel('Frequency')
    ylabel('Power Spectrum Magnitude (dB)')
end