www.gusucode.com > 生物信号工具箱 > 生物信号工具箱/生物信号工具箱/Biosignal/Biosignal/Graphs/fftplot.m
function [varargout]=fftplot(varargin) % %[Pxx,f]=fftplot(x,nfft,Fs,scale,show,clr) % % % _______________________________________________________________________________ % Copyleft (l) 2010 by Ikaro Silva, All Rights Reserved % Contact Ikaro Silva (ikarosilva@ieee.org) % % This library (Biosignal Toolbox) is free software; you can redistribute % it and/or modify it under the terms of the GNU General Public License as published by % the Free Software Foundation; either version 2 of the License, or % (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program; if not, write to the Free Software % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA % 02111-1307 USA % % _______________________________________________________________________________ %Initialize Input Parameter Names and Default values tmp=2^(nextpow2(varargin{1})); %Use next highest power of 2 greater than or equal to length(x) to calculate FFT. in_params={'x,nfft,Fs,scale,show,clr'}; in_params_default={varargin{1},tmp,[],1,[],[0.5 0.5 0.5]}; get_varargin(in_params,in_params_default,varargin); % Take fft, padding with zeros so that length(fftx) is equal to nfft fftx = fft(x,nfft); % Calculate the numberof unique points NumUniquePts = ceil((nfft+1)/2); % FFT is symmetric, throw away second half fftx = fftx(1:NumUniquePts); % Take the magnitude of fft of x and scale the fft so that it is not a function of % the length of x mx = abs(fftx)/length(x); % Take the square of the magnitude of fft of x. mx = mx.^2; % Since we dropped half the FFT, we multiply mx by 2 to keep the same energy. % The DC component and Nyquist component, if it exists, are unique and should not % be mulitplied by 2. if rem(nfft, 2) % odd nfft excludes Nyquist point mx(2:end) = mx(2:end)*2; else mx(2:end -1) = mx(2:end -1)*2; end if(scale) %normalize the power to peak at 1 mx=mx./max(mx); end % This is an evenly spaced frequency vector with NumUniquePts points. f = linspace(0,1,NumUniquePts); %normalized frequency %Convet frequency vector to Hertz if Fs is passed in if(~isempty(Fs)) f=f.*(Fs/2); end %Convert power to dB mx=10*log10(mx); % Generate the plot, title and labels. if(~isempty(show)) switch lower(show) case 'logylogx' semilogx(f,mx,'Color',clr); case {'logylinx','logxliny'} plot(f,mx,'Color',clr); case 'linylinx' plot(f,10.^(mx./10),'Color',clr); case {'linylogx','linxlogy'} semilogx(f,10.^(mx./10),'Color',clr); end title('Power Spectrum (dB)'); ylabel('Power (dB)'); grid on if(isempty(Fs)) xlabel('Normalized Frequency'); else xlabel('Frequency'); end end %Pass back the outputs if requested if(nargout >= 1) varargout(1)={mx}; if(nargout >= 2) varargout(2)={f}; end end