Welch and MEM power spectrum estimation (Signal processing) functions
Peter Vernon Lanspeary <[email protected]> Thu, 15 Jun 2006 20:04:41 +0930
| Newsgroups | gmane.comp.gnu.octave.sources |
|---|---|
| Message-ID | <[email protected]> |
--7AUc2qLy4jB3hD7Z
Content-Type: text/plain; charset=us-ascii
Content-Disposition: inline
Here are some (attached) functions for estimating power spectra of vector data.
1) welch_psd.m -- (Welch 1967) averaged-periodogram method of spectral
estimation: similar to Matlab's pwelch
2) burg_filter.m -- constructs autoregressive generating filter of the data
by the Burg (1968) lattice-filter algorithm
3) mem_psd.m -- maximum entropy method of spectral estimation (Burg 1968
et seq); uses filter coefficients from burg_filter.m
mem_psd.m and burg_filter.m perform a similar function to
Matlab's pburg
These have been slightly altered since testing. Please inform me of bugs.
Thanks
Peter Lanspeary
--7AUc2qLy4jB3hD7Z
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename="burg_filter.m"
function [coeffs,residual] = burg_filter(data,poles,stop_crit)
%%
%% [coeffs,residual] = burg_filter(data,poles,stop_crit)
%%
%% Please report bugs to [email protected]
%%
%% Calculates coefficients of a moving average (MA) whitening filter
%% using the lattice-filter of Burg (1968) and also described by Kay
%% and Marple (1981). This filter reduces the data to white noise.
%% Two other filters can be built from this whitening filter:
%% (1) a moving average linear prediction filter -- by removing the
%% first (zero-lag) coefficient and negating the others.
%% (2) an autoregressive generating (AR) filter -- by using "residual"
%% as the sole MA coefficient and using "coeffs" to supply the AR
%% coefficients of an ARMA filter.
%% The power spectrum of the ARMA filter is an estimate of the maximum
%% entropy power spectrum of the data.
%%
%% Arguments:
%% data %% [real vector] sampled data
%% poles %% [integer scalar] required number of poles of AR filter
%% stop_crit %% [string] if 'FPE' or 'AIC', the FPE or AIC criterion
%% %% (Kay & Marple, 1981) are used to override the 'poles'
%% %% argument so the filter does not grow too long.
%%
%% Returned values:
%% coeffs %% [real vector] list of M=(poles+1) moving-average filter
%% %% coefficients; for data input x(n) and
%% %% white noise output e(n), the filter is
%% %% M
%% %% e(n) = SUM coeffs(k).x(n-k)
%% %% k=0
%% %% N.B. the first element of "coeffs" is the zero-lag
%% %% coefficient, which always has a value coeffs(1)=1.
%%
%% residual %% mean square of residual (white) noise from filter
%%
%% REFERENCES
%% John Parker Burg (1968):
%% "A new analysis technique for time series data",
%% NATO advanced study Institute on Signal Processing with Emphasis on
%% Underwater Acoustics, Enschede, Netherlands, Aug. 12-23, 1968.
%%
%% Steven M. Kay and Stanley Lawrence Marple Jr. (1981):
%% "Spectrum analysis -- a modern perspective",
%% Proceedings of the IEEE, Vol 69, pp 1380-1419, Nov., 1981
%%
1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Copyright 2006 (C) Peter V. Lanspeary
%%
%% burg_filter.m 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, 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
%% (http://www.gnu.org/copyleft/gpl.html) for more details.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%
%% sanity checks
if ( nargin < 2 )
error( 'burg_filter(data,poles): Need at 2 args.' );
elseif ( ~ isvector(data) || ~ isreal(data) || length(data) < 3 )
error( 'burg_filter:error: arg 1 (data) must be real vector of length >3.' );
elseif ( ~isscalar(poles) || ~isreal(poles) || fix(poles)~=poles || poles<=0.5)
error( 'burg_filter:error: arg 2 (poles) must be positive integer.' );
elseif ( floor(poles+0.5) > length(data)-2 )
error( 'burg_filter:error: arg 2 (poles)+2 must be less than data length' );
elseif ( nargin>2 && ( ~ischar(stop_crit) || size(stop_crit,1)>1 ) )
error( 'burg_filter:error: arg 3 (stop_crit) must be one string' );
else
%% end of sanity checks
%%
coeffs = zeros(1,poles);
%%
%% Storage of forward and backward prediction errors is a little tricky.
%% Because the forward error e(i) is always combined with the lagged
%% backward error b(i-1), e(1) and b(n) are never used, and therefore are
%% never stored. Not storing unused data makes the calculation of the
%% reflection coefficient look much cleaner :)
data = reshape(data,1,[]);
N = size(data, 2);
forw_err = data(2:N);
back_err = data(1:N-1);
residual = sumsq(data)/N;
%%
%% Decide if FPE or AIC criteria will be applied to stop before
%% getting to the specified number of filter poles.
if ( nargin > 2 )
use_FPE = strcmp(stop_crit,'FPE');
use_AIC = strcmp(stop_crit,'AIC');
else
use_FPE = 0;
use_AIC = 0;
endif
new_criterion = residual;
old_criterion = 2 * new_criterion;
for k = 1:poles
%%
%% reflection_coeff = -2 * E(e(i)*b(i-1)) / ( E(e(i)^2) + E(b(i-1)^2) )
refl_coeff= -2 *forw_err * back_err' / (sumsq(forw_err) + sumsq(back_err));
%% Levinson-Durbin recursion for residual
new_residual = residual * ( 1.0 - refl_coeff^2 );
if ( k > 1 )
%%
%% Apply the FPE or AIC criterion and stop if the FPE or AIC is
%% increasing rather than decreasing.
%% Do it before we update the old filter "coeffs" and "residual".
if ( use_FPE )
old_criterion = new_criterion;
new_criterion = new_residual * ( N + k + 1 ) / ( N - k - 1 );
if ( new_criterion > old_criterion )
break;
endif
elseif ( use_AIC )
old_criterion = new_criterion;
new_criterion = log(new_residual) + 2 * ( k + 1 ) / N;
if ( new_criterion > old_criterion )
break;
endif
endif
%% Update filter "coeffs" and "residual".
%% Use Levinson-Durbin recursion formula.
coeffs = [ prev_coeffs+refl_coeff.*prev_coeffs(k-1:-1:1), refl_coeff ];
residual = new_residual;
else
coeffs = refl_coeff;
residual = new_residual;
end
if ( k < poles )
prev_coeffs = coeffs;
%% calculate new prediction errors (by recursion):
%% e(p,i) = e(p-1,i) + k * b(p-1,i-1) i=2,3,...n
%% b(p,i) = b(p-1,i-1) + k * e(p-1,i) i=2,3,...n
%% remember e(p,1) is not stored, so don't calculate it; make e(p,2)
%% the first element in forw_err. b(p,n) isn't calculated either.
nn = length(forw_err);
forw_new = forw_err(2:nn) + refl_coeff .* back_err(2:nn);
back_err = back_err(1:nn-1) + refl_coeff .* forw_err(1:nn-1);
forw_err = forw_new;
end
end
%% end of for loop
coeffs = [1 coeffs];
end
end
--7AUc2qLy4jB3hD7Z
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename="mem_psd.m"
function [ret_psd,ret_freq]=mem_psd(ar_coeffs,residual,freq,sample_f,method)
%%
%% [ret_psd,ret_freq] = mem_spect(ar_coeffs,residual,freq,sample_f,method)
%%
%% Please report bugs to [email protected]
%%
%% Calculates the power spectrum of the autoregressive filter
%%
%% M
%% x(n) = residual.e(n) + SUM ar_coeffs(k).x(n-k)
%% k=1
%% where x(n) is the filter output and e(n) is white noise input.
%% If the "freq" argument is a vector (of frequencies) the spectrum is
%% calculated using the polynomial method and the "method" argument is
%% ignored. For scalar "freq", an integer power of 2, or "method='FFT'",
%% causes the spectrum to be calculated by FFT. Otherwise, the spectrum
%% is calculated as a polynomial. Note that it is more computationally
%% efficient to use the FFT method if length of the filter is not much
%% smaller than the number of frequency values. The spectrum is scaled so
%% that spectral energy between zero frequency and the Nyquist frequency
%% is the same as the time-domain energy (i.e. mean square of the signal).
%%
%% Arguments:
%% ar_coeffs %% [real vector] list of M=(order+1) autoregressive filter
%% %% N.B. the first element of "coeffs" is the zero-lag
%% %% coefficient, which always has a value coeffs(1)=1.
%%
%% residual %% [real scalar] the moving-average coefficient of the AR
%% %% filter.
%% freq %% [real vector] frequencies at which power spectral density
%% %% is calculated
%% %% [integer scalar] number of frequency values (uniformly
%% %% distributed from zero to the Nyquist frequency) at
%% %% which spectral density is calculated. [default=256]
%% sample_f %% [real scalar] sampling frequency (Hertz) [default=1]
%% method %% [string] controls the method of evaluation if "freq"
%% %% is a scalar ---
%% %% method="FFT": use FFT to calculate power spectrum.
%% %% method="poly": use polynomial method
%% %% [default="poly" unless number of frequencies is an
%% %% integer power of 2
%% freq, sample_f and method are optional.
%% freq and sample_f may be empty.
%%
%% Returned values:
%% If return values are not required by the caller, the spectrum
%% is plotted and nothing is returned.
%% ret_psd %% [real vector] power-spectrum estimate
%% ret_freq %% [real vector] frequency values
%%
%% REFERENCES
%% William H. Press and Saul A. Teukolsky and William T. Vetterling and
%% Brian P. Flannery",
%% "Numerical recipes in C, The art of scientific computing", 2nd edition,
%% Cambridge University Press, 2002 --- Section 13.7.
%% N.B. The algorithm in Press et al. expects prediction-filter coefficients.
%% mem_spect requires the whitening-filter coefficients.
%%
1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Copyright 2006 (C) Peter V. Lanspeary
%%
%% mem_spect.m 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, 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
%% (http://www.gnu.org/copyleft/gpl.html) for more details.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% sanity checks
if ( nargin > 2 )
user_freqs = isvector(freq) && length(freq)>1;
end
if ( nargin < 2 )
error( 'mem_spect(ar_coeffs,residual,freq,sample_f,method) needs >=2 args');
elseif ( ~ isvector(ar_coeffs) || ~ isreal(ar_coeffs) || length(ar_coeffs)<2 )
error( 'mem_spect: error: arg 1 (ar_coeffs) not real vector of length>=2.' );
elseif ( ~ isscalar(residual) || ~ isreal(residual) || residual <= 0 )
error( 'mem_spect: error: arg 2 (residual) is not real scalar > 0' );
elseif ( nargin > 2 && ~isscalar(freq) && ~user_freqs && ~isempty(freq) )
error( 'mem_spect: error: arg 3 (freq) is not vector or scalar.');
elseif ( nargin > 2 && ~user_freqs && ~isempty(freq) && ...
( ~isreal(freq) || fix(freq)~=freq || freq <= 2 || freq >= 1048576 ) )
error( 'mem_spect: error: arg 3, (freq) is not integer >=2, <=1048576.)' );
elseif ( nargin > 2 && user_freqs && ~isempty(freq) && ...
( ~all(isreal(freq)) || any(freq<0) ) )
error( 'mem_spect: error: arg 3, (freq) values must be real and >=0' );
elseif ( nargin > 3 && ~isempty(sample_f) && ...
( ~isscalar(sample_f) || ~isreal(sample_f) || sample_f<0 ) )
error( 'mem_spect: error: arg 4, sample_f must be real scalar >0.' );
elseif ( nargin > 4 && ~ ischar(method) && ...
~ strcmp(method,'FFT') && ~ strcmp(method,'poly') )
error( 'mem_spect: error: arg 5, (method) must be "FFT" or "poly".' );
else
%% end of sanity checks
%%
%% define the frequencies
if ( nargin < 3 || ( nargin >= 4 && isempty(freq) ) )
freq = 256;
user_freqs = 0;
end
if ( nargin < 4 || ( nargin >= 5 && isempty(sample_f) ) )
sample_f = 1.0;
end
if ( user_freqs )
len_freq = length(freq);
else
len_freq = freq;
freq = [0:len_freq] * sample_f / 2 / len_freq;
end
%%
%% decide which method to use
is_power_of_2 = rem(log(len_freq),log(2))<10.*eps;
force_FFT = nargin>4 && strcmp(method,'FFT');
force_poly = nargin>4 && strcmp(method,'poly');
use_FFT = ~user_freqs && ( ( ~ force_poly && is_power_of_2 ) || force_FFT);
%%
%% and do it
len_coeffs = length(ar_coeffs);
if ( use_FFT )
fft_out = fft( [ ar_coeffs zeros(1,len_freq*2-len_coeffs ) ] );
fft_out = fft_out(1:len_freq+1);
else %% do not use FFT
two_pi_i = (0+i) * 2 * pi / sample_f;
fft_out = ar_coeffs * exp( two_pi_i * [1:len_coeffs]' * freq );
end
spectrum = 2.0 * residual / sample_f ./ ( fft_out .* conj(fft_out) );
%%
if ( nargout <= 0 )
loglog(freq,spectrum);
else
ret_psd = spectrum;
ret_freq = freq;
end
end
end
--7AUc2qLy4jB3hD7Z
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename="welch_psd.m"
function ...
[ret_psd,ret_freq]=welch_psd(signal,window,sample_f,overlap,padding,qual)
%%
%% [ret_psd,ret_freq]=welch_psd(signal,window,sample_f,overlap,padding,qual)
%%
%% Estimate power spectrum of a time-series signal by the periodogram
%% (FFT) method. The "signal" is divided into segments of length
%% equal to the length of "window", and each segment is multiplied by
%% "window" before (optional) zero-padding and calculating its FFT.
%% The power spectrum is the mean square of the FFTs, scaled so that area
%% under the power spectrum is the same as the mean square of the signal.
%%
%% Arguments
%% signal %% [real vector] time-series signal
%%
%% window %% [real vector] of window-function values between 0 and 1; the
%% %% signal segment has the same length as the window, or
%% %% [integer scalar] length of each signal segment (and
%% %% unpadded FFT);
%% %% default value is sqrt(length(x)) rounded down to the
%% %% nearest integer power of 2, and forces qual="sloppy"
%%
%% sample_f %% [real scalar] sampling frequency (Hertz); default=1.0
%%
%% overlap %% [real scalar] segment overlap factor
%% %% 0 <= overlap < 1, default is zero, 0.5 is "industry standard"
%%
%% padding %% [integer scalar] number of samples of zero padding (per FFT)
%% %% default is zero
%%
%% qual %% [string] Quality specifier
%% %% "sloppy" FFT length is rounded up to the nearest integer
%% %% power of 2, segment length is rounded up unless
%% %% the "window" arg is specified as a vector
%% %% FFT length is adjusted after addition of padding
%% %% default is to use exactly the segment lengths and padding
%% %% lengths (hence FFT length) specified in argument list
%% All but the first argument are optional.
%% Any but the first and last may be empty.
%%
%% Returned values:
%% If return values are not required by the caller, the spectrum
%% is plotted and nothing is returned.
%% ret_psd %% [real vector] power-spectrum estimate
%% ret_freq%% [real vector] frequency values
%% %% the length of both returned vectors is [length_of_FFT]/2+1
%%
%% REFERENCE
%% Peter D. Welch (June 1967):
%% The use of fast Fourier transform for the estimation of power spectra:
%% a method based on time averaging over short, modified periodograms.
%% IEEE Transactions on Audio Electroacoustics, VOl AU-15(6), pp 70-73
%%
1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Copyright 2006 (C) Peter V. Lanspeary
%%
%% welch_psd.m 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 Soft-
%% ware Foundation; either version 2, 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
%% (http://www.gnu.org/copyleft/gpl.html) for more details.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
if ( nargin >0 )
x_len = length(signal);
is_win = 0;
if ( nargin > 1 )
if ( isscalar(window) )
is_win = 1;
elseif ( isvector(window) )
is_win = length(window);
end
end
end
%%
%% SANITY CHECKS
if ( nargin <= 0 )
error( 'welch_psd(signal,window,overlap,padding): Need at least 1 arg.' );
elseif ( ~isvector(signal) || ~ isreal(signal) )
error( 'welch_psd: error: arg 1 (signal) must be real vector.' );
elseif ( nargin > 1 && ~isempty(window) && ~is_win )
error( 'welch_psd: error: arg 2 must be window vector or segment length.' );
elseif ( nargin > 1 && is_win==1 && ( ~isreal(window) || ...
fix(window)~=window || window<4 || x_len<window ) )
error( 'welch_psd: error: arg 2, window not integer, >4 & <=length(data).' );
elseif ( nargin > 1 && is_win>1 && ( ~all(isreal(window)) || any(window<0) ) )
error( 'welch_psd: error: arg 2, window values must be real and >=0' );
elseif ( nargin > 2 && ~isempty(sample_f) && ...
( ~isscalar(sample_f) || ~isreal(sample_f) || sample_f<0 ) )
error( 'welch_psd: error: arg 3, sample_f must be real scalar >0.' );
elseif ( nargin > 3 && ~isempty(overlap) && ...
( ~isscalar(overlap) || ~isreal(overlap) || overlap<0 || overlap>=1 ) )
error( 'welch_psd: error: arg 4, overlap not between 0 and 1' );
elseif ( nargin > 4 && ~isempty(padding) && ...
( ~isscalar(padding) || ~isreal(padding) || ...
fix(padding)~=padding || padding<0 ) )
error( 'welch_psd: error: arg 5, padding not integer >=0' );
elseif ( nargin > 5 && ( ~ischar(qual) || size(qual,1)>1 ) )
error( 'welch_psd: error: arg 6, qual must be one string' );
else
%% end of sanity checks
%%
log_two = log(2);
nearly_one = 0.99999999999;
is_sloppy = ~is_win || ( is_win==1 && nargin>5 && strcmp(qual,'sloppy') );
%%
%% calculate/adjust segment length
if ( ~is_win )
seg_len = 2 ^ ceil( log( sqrt(x_len) ) * nearly_one / log_two );
window = 1;
win_meansq = 1;
elseif ( is_win==1 )
seg_len = window;
window = 1;
win_meansq = 1;
else
window = reshape(window,1,[]);
seg_len = length(window);
win_meansq = sumsq(window) / seg_len;
end
%%
if ( nargin<3 || ( nargin>=3 && isempty(sample_f) ) )
sample_f = 1;
end
if ( nargin<4 || ( nargin>=4 && isempty(overlap) ) )
overlap = 0;
end
if ( nargin<5 || ( nargin>=5 && isempty(padding) ) )
padding = 0;
end
if ( x_len < seg_len )
error( 'welch_psd: error: Signal is shorter than segment/window' );
end
%%
%% calculate FFT length
fft_len = seg_len + padding;
if ( is_sloppy )
fft_len = 2 ^ ceil( log( fft_len ) * nearly_one / log_two );
end
%%
%% AVERAGE THE SIGNAL
%%
%% Use the same signal data as the periodogram to take into account the
%% overlap and any unused data.
%% N.B. that removing the mean for ALL the data makes the spectral energy
%% a bit too high, removing a "local/segment" mean value from each segment
%% makes spectral energy too low. What is the CORRECT method? IDNK.
%%
overlap = fix(seg_len * overlap);
avg_signal = 0;
n_fft = 0;
for start_seg = [1:seg_len-overlap:x_len-seg_len+1]
avg_signal = avg_signal + sum(signal(start_seg:start_seg+seg_len-1));
n_fft = n_fft +1;
end
avg_signal = avg_signal / n_fft / seg_len;
%%
%% CALCULATE PERIODOGRAM
%%
%% When computing the FFT of real data, it is possible to increase the
%% speed of the algorithm by
%% x = signal(start_seg:2:end_seg) + i * signal(start_seg+1:2:end_seg);
%% and then unscrambling the vector returned by fft(). The length of
%% the decreaes by a factor of 2. However, it is quite possible that
%% fft() recognises real input and uses this (or a similar) trick anyway.
%%
psd_len = floor(fft_len/2)+1;
psd = zeros(1,psd_len);
x = zeros(1,fft_len);
n_fft = 0;
for start_seg = [1:seg_len-overlap:x_len-seg_len+1]
x(1:seg_len) = window .* \
(signal(start_seg:start_seg+seg_len-1) - avg_signal );
X = fft(x);
n_fft = n_fft +1;
psd = psd + abs( X(1:psd_len) ).^2;
end
psd = psd / ( n_fft * seg_len * sample_f * win_meansq / 2 );
freq = [0:psd_len-1] * ( sample_f / fft_len );
if ( nargout <= 0 )
loglog(freq,psd);
else
ret_psd = psd;
ret_freq = freq;
end
end
end
--7AUc2qLy4jB3hD7Z
Content-Type: text/plain; charset="us-ascii"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Content-Disposition: inline
_______________________________________________
Octave-sources mailing list
[email protected]
https://www.cae.wisc.edu/mailman/listinfo/octave-sources
--7AUc2qLy4jB3hD7Z--