lpsvd

Linear Prediction and Singular Value Decomposition

Syntax

predictedSpectrum = lpsvd(spectrum, time, method, order)
[predictedSpectrum, PredictionParameters] = lpsvd(...)

Description

lpsvd estimates the parameters for a linear combination of exponentially damped sinusoids needed to fit the input spectrum through linear prediction singular value decomposition or related methods.

The model for the exponentially damped sinusoid is:

y(time) = amp * exp(1i*phase) * exp(time * (1i*2*pi*freq - damp) )

The function requires 2-4 input parameters: the spectrum to be fit and its corresponding time vector are necessary; the method and order if not provided will revert to the defaults.

The third input determines the method used to fit the model to the spectrum, the three possible inputs include:

'kt': Kumaresan and Tufts linear prediction and singular value decomposition algorithm
'ss': State space and singular value decomposition method (default)
'tls': Hankel total least squares method

The fourth input is an initial estimate for the number of damped sinusoids present in the spectrum. If an integer is provided that will be used otherwise the order can be estimated via:

'mdl': minimum description length (default)
'aic': Akaike information protocol

However these two methods are known to underestimate the number of components.

Examples

Generate an example set of damped sinusoids and add noise

clear
dt = 4/1000; % µs
n = 1024;
time = (0:n-1)*dt;

T2 = 0.4;  % µs
w = -50:10:50;  % MHz
A = -1.1:.2:0.9;

y = exp(1i*0)*exp(-time/T2).*(A*exp(1i*2*pi*w'*time));
yn = complex(real(y)+ 0.3*randn(1,n),imag(y) + 0.3*randn(1,n));

yf = fftshift(fft(y));
ynf = fftshift(fft(yn));
xf = fdaxis(dt,n);

subplot(2,1,1)
plot(time,imag(y),time,imag(yn)); xlabel('time (ns)')
subplot(2,1,2)
plot(xf,real(yf),xf,real(ynf)); xlabel('frequency (MHz)'); xlim([-100 100])

Perform lpsvd on the noisy data

yp = lpsvd(yn,time,'ss','aic');
ypf = fftshift(fft(yp));

subplot(2,1,1)
plot(time,imag(yn),time,imag(yp)); xlabel('time (ns)')
subplot(2,1,2)
plot(xf,real(ynf),xf,real(ypf)); xlabel('frequency (MHz)'); xlim([-100 100])

References

The methods implemented in lpsvd are described in