ezfft.m

Computes/displays the power spectrum of a real time series

function [w,e] = ezfft(varargin)
%EZFFT  Easy to use Power Spectrum
%   EZFFT(T,U) plots the power spectrum of the U(T) time series, where T
%   is the 'time' and U is a real signal.  If T is a scalar, then it is
%   interpreted as the sampling time of the signal U.  If T is a vector,
%   then it is interpreted as the time itself. In this latter case, T must
%   be equally spaced (as obtained by LINSPACE for instance), and it must
%   have the same length as U. If T is not specified, then a sampling time
%   of unity (1 second for instance) is taken.
%
%   [W,E] = EZFFT(T,U) returns the power spectrum E(W), where E is the
%   energy density and W the pulsation 'omega'.  W is *NOT* the frequency.
%   The frequency is W/(2*pi).
%
%   The length of the vectors W and E is N/2, where N is the length of U
%   (this is because U is assumed to be a real signal.) If N is odd, the
%   last point of U (and T) is ignored.
%
%       W(1) is always 0.  E(1) is the energy density of the average of U.
%       W(2) is the increment of pulsation, Delta W, given by 2*PI/Tmax
%       W(end), the highest measurable pulsation, is PI/DT, where DT is the
%          sampling time (Nyquist theorem).
%
%   Energy conservation (Parseval theorem)
%   For every signal U, one has the following property:
%      MEAN(U.^2) == SUM(E)*W(2)
%   (this works only if no appodization is applied. In case of an
%   appodization, some energy is lost, so the spectral energy is lower than
%   the actual one.)
%
%   EZFFT(..., 'Property1', 'Property2', ...) specifies the properties:
%    'hann'    applies a Hann appodization window to the data (reduces
%              aliasing).
%    'disp'    displays the spectrum (by default if no output argument)
%    'freq'    the frequency f is displayed instead of the pulsation omega
%              (this applies for the display only: the output argument
%              remains the pulsation omega, not the frequency f).
%    'handle'  returns a handle H instead of [W,E] - it works only if the
%              properties 'disp' is also specified. The handle H is useful      
%              to change the line properties (color, thickness) of the plot.
%
%   Example 1:
%      t = linspace(0,400,2000);
%      u = 0.2 + 0.7*sin(2*pi*t/47) + cos(2*pi*t/11);
%      [w,e] = ezfft(t,u,'disp');
%
%   Example 2:
%      h = ezfft(t,u,'disp','handle');
%      set(h,'Color','red');
%
%   Example 3:
%      [w,e] = ezfft(t,u,'hann');
%      loglog(w,e,'b*');
%
%   See also FFT

%   F. Moisy, moisy_at_fast.u-psud.fr
%   Revision: 1.00,  Date: 2008/11/19


error(nargchk(1,inf,nargin));

if nargin==1
    u=varargin{1};
    dt=1;
else
    dt=varargin{1};
    u=varargin{2};
end

if length(dt)>1
    if length(u)~=length(dt)
        error('Vectors T and U should be of equal length.');
    end
    dt = abs(dt(2)-dt(1));
end

% works only with vectors of even length
if mod(length(u),2)~=0;
    u=u(1:(end-1));
end

u=real(u);

% works with line vector (otherwise transpose it)
if size(u,2)==1
    u=u';
end

% signal appodization with a Hann window
if any(strcmpi(varargin,'hann'))
    u=u.*hann(length(u))';
end

n = length(u)/2;  % length of the output fft
w = linspace(0,pi,n)/dt;   % output pulsation

e = abs(fft(u)).^2;
e = 2*e(1:n) / (2*n)^2 / w(2);  % normalisation
e(1)=e(1)/2;    % mode w=0 (twice energy is the previous line)

% display
if nargout==0 || any(strncmpi(varargin,'disp',1))
    if any(strncmpi(varargin,'freq',1))
        hh=loglog(w/(2*pi),2*pi*e);
        xlabel('f');
        ylabel('E(f)');
    else
        hh=loglog(w,e);
        xlabel('\omega');
        ylabel('E(\omega)');
    end
end

if nargout>0 && any(strncmpi(varargin,'handle',4)) && any(strncmpi(varargin,'disp',1))
    w=hh;
    clear e
    return
end

if nargout==0
    clear w e
end


function y=hann(n)
%HANN  Window Hann function
%   Y = HANN(N) returns a N-point Hann function (usually N is a power of
%   2), ie a discretization of Y(X) = (1 - COS (2*PI*X))/2.
%
%   The HANN function is used to window (apodize) finite samples of signal
%   in order to avoid high frequency oscillations in Fourier transform.
%
%   Example:  plot(hann(128));
%
%   F. Moisy
%   Revision: 1.03,  Date: 2006/03/03
%
%   See also FFT.

% History:
% 2004/09/17: v1.00, first version.
% 2005/23/02: v1.01, cosmetics
% 2005/09/03: v1.02, id.
% 2006/03/03: v1.03, bug fixed 0:n-1

error(nargchk(1,1,nargin));
y=0.5*(1-cos(2*pi*(0:n-1)/(n-1)))';