ohspeci.m

一种间接的估计信号的1.5维谱的方法

function [OHspec,waxis] = ohspeci(y,nlag,nsamp,overlap,nfft, wind)
% OHSPECI one and half spectrum estimation using the indirect method. 
%	[OHspec,waxis] = ohspeci (y,nlag)
%	[OHspec,waxis] = ohspeci (y,nlag,nsamp,overlap,nfft, wind)
%	y       - data vector or time-series   
%	nlag    - number of lags to compute [must be specified] 
%	nsamp   - samples per segment    [default: length of y]
%	overlap - percentage overlap     [default = 0] 
%	nfft    - FFT length to use      [default = 128] 
%	wind    - window function to apply: 
%	      if wind=0, the hanning window is applied (default); 
%	      otherwise the hamming window with unity values is applied. 
%	OHspec  - estimated one and half spectrum
%	waxis   - frequency-domain axis associated with the one and half spectrum.  


% ------------------- parameter checks ---------------------

    ly = length (y); 
    if (exist('overlap') ~= 1)   overlap = 0;          end 
    overlap = min(99, max(overlap,0)); 
    if (exist('nsamp') ~= 1)     nsamp   = ly;         end
    if (nsamp > ly | nsamp <= 0) nsamp   = ly;         end
    if (exist('nfft') ~= 1)      nfft    = 128;        end
    if (nfft <= 0)               nfft    = 128;        end 
    if (exist('wind') ~= 1)      wind    = 0;          end

    nlag = min(nlag, nsamp-1); 
    if (nfft  < 2*nlag+1)   nfft = 2^nextpow2(2*nlag+1);  end 
    
% ---------------- create the lag window --------------------

    if (wind == 0) 
         window = hann(2*nlag+1);
     else
         window = hamming(2*nlag+1); 
     end 

% ------------- one and half spectrum estimation -------------
    
    overlap  = fix(nsamp * overlap / 100); 
    nadvance = nsamp - overlap; 
    nrecord  = fix ( (ly - overlap) / nadvance );
    
    OHspec = zeros(1,nfft); 
    C = zeros(1,2*nlag+1);
    ind = [1:nsamp];
    for k = 1:nrecord
        x = y(ind);
        x = x - mean(x);
        ind = ind + nadvance; 
        c = zeros(1,2*nlag+1);
        for i = -nlag : nlag
            M1 = max(0,-i);
            M2 = min(nsamp-1,nsamp-1-i);
            r = 0;
            for j = M1:M2 
                r = r + x(j+1) * x(i+j+1) * x(i+j+1);
            end;
            c(i+nlag+1) = r/nsamp;
        end;
        C = C+c;
    end;
    C = C/nrecord;
    C = C.*window';
    fft_C= fft(C,nfft);
    fft_C = fftshift(fft_C);
    OHspec = abs(fft_C);
    OHspec = OHspec/max(OHspec);

% ---------------- create the waxis --------------------

    if (rem(nfft,2) == 0) 
        waxis = [-nfft/2:(nfft/2-1)]/nfft; 
    else
        waxis = [-(nfft-1)/2:(nfft-1)/2]/nfft; 
    end ;  

% % ---------------- plot OHspec ----------------  
%     figure;
%     plot(waxis,OHspec);
%     title('One and half spectrum estimation')
%     xlabel('Normalized Frequency');
%     ylabel('Normalized Amplitude') 
% % ----------------------------------------------
 return;
% ----------- End od function file --------------