Commit 3dfad0f4 authored by Philip Roy Charlton's avatar Philip Roy Charlton

Removed outdated median PSD files

parent f57b5d1d
clear all;
nFFT = 64;
nOverlap = nFFT/2;
%nOverlap = 0;
windowFn = ones(nFFT, 1);
Fs = 2048.0;
% How many segments
M = 1024;
nx = nFFT + (M-1)*(nFFT - nOverlap);
x = rand(nx, 1);
tic;
[pmed, fmed] = medianPSD(x, windowFn, nOverlap, nFFT, Fs, 'median');
%[pmed, fmed] = medianPSD(x, windowFn, nOverlap, nFFT, Fs, 'mean');
pmed = pmed.';
fmed = fmed.';
toc;
tic;
%[pavg, favg] = pwelch(x, windowFn, nOverlap, nFFT, Fs, 'twosided');
[pavg, favg] = pwelch(x, windowFn, nOverlap, nFFT, Fs);
pavg = pavg.';
favg = favg.';
toc;
pmed
pavg
\ No newline at end of file
%
% medianPSD --- estimate the PSD of a time-series using the median method
%
% [p, f] = medianPSD(x, windowFn, nOverlap, nFFT, Fs, method)
%
% estimates the power spectral density of x by calculating periodograms
% of overlapping segments in a times-series x and taking the median value
% in each frequency bin, adjusted by a bias factor. The parameters play the
% same roles as they do in Welch's method implemented in pwelch(), however
% in pwelch() the periodograms are averaged.
%
% The sample median is being used to estimate for the population mean mu of
% the exponential distribution in each frequency bin becuase it is more robust
% to outliers in the PSD of real data. The relationship is
%
% mu = (population median)/ln(2)
%
% however, when using the sample median to estimate mu, the bias factor depends on
% sample size, approaching 1/ln(2) for large samples. This function incorporates
% the bias correction factor calculated in Allen et. al.
%
% https://arxiv.org/abs/gr-qc/0509116
%
% Another source of bias comes from the fact that the overlapping segments are not
% independent. This can be mitigated by only taking the median of two sets of
% non-overlapping segments and then taking the average of the two PSDs, as described
% in Allen et. al., but that has not been implemented in this function.
%
% Input:
%
% x - a time-series
% windfowFn - a window function, which must be a vector of length nFFT
% nOverlap - the number of samples that will overlap in each segment of x
% nFFT - the FFT length used for periodograms
% Fs - the sampling frequency (Hz)
% (optional) method - a string, either 'median' or 'mean', to determine
% which method will be used. For 'mean', the results should be identical
% to those given by pwelch(). The default is 'median'.
%
% Output:
%
% p - the PSD in [x]^2/Hz
% f - bin centre frequencies (Hz)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [p, f]=medianPSD(x, windowFn, nOverlap, nFFT, Fs, method)
% Set default values
if (nargin < 6)
method = 'median';
end;
% If there are M segments then the total length is
%
% nx = nFFT + (M-1)*(nFFT - nOverlap)
%
% so for length nx we need
%
% M = (nx - nFFT)/(nFFT - nOverlap) + 1
%
% Check arguments
nx = length(x);
wn = length(windowFn);
if (wn ~= nFFT)
error('Window length is different to FFT length');
end
if (nOverlap < 0)
error('Overlap is negative');
elseif (nOverlap >= nFFT)
error('Overlap is >= FFT length');
end;
M = ceil((nx - nFFT)/(nFFT - nOverlap)) + 1;
wt = norm(windowFn).^2;
dataIsReal = isreal(x);
% Weighting factor for the PSD
if (dataIsReal)
fac = 2/(wt*Fs);
else
fac = 1/(wt*Fs);
end;
data = zeros(nFFT, M);
for k = 1:M
idx = (k - 1)*(nFFT - nOverlap) + 1;
data(:, k) = windowFn.*x(idx:idx+nFFT-1);
end;
data = fft(data, [], 1);
if (dataIsReal)
data = abs(data(1:nFFT/2+1, :)).^2;
else
data = abs(data).^2;
end;
if (strcmp(method, 'median'))
p = (fac/medianBias(M))*median(data, 2);
elseif (strcmp(method, 'mean'))
p = fac*mean(data, 2);
else
error('Unknown PSD method');
end;
if (dataIsReal)
% One-sided psd
f = (Fs/nFFT)*[0:nFFT/2].';
p(1) = 0.5*p(1);
p(end) = 0.5*p(end);
else
% Two-sided psd
f = (Fs/nFFT)*[0:nFFT-1].';
end;
return;
function biasFac=medianBias(M)
L = [1:M];
biasFac = sum((-1).^(L+1).*(1./L));
return;
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment