View Raw SPL
/*****************************************************************************
* *
* INVPSD.SPL Copyright (C) 2006-2025 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Constructs a time series given a PSD *
* *
* Revisions: 17 May 2006 RRR Creation *
* 3 Mar 2025 RRR frequency domain interpolation *
* *
*****************************************************************************/
#if @HELP_INVPSD
INVPSD
Purpose: Constructs a time series with a PSD that matches the PSD input.
Syntax: INVPSD(psd, N, rpflag, "method")
psd - A series, the input PSD to match.
N - Optional. An integer, the time series output length.
Defaults to 2 * (length(psd) - 1).
rpflag - Optional. An integer.
0: Use zero phase to construct time series.
1: Use random phase to construct time series (default).
"method" - Optional. A string, the frequency domain interpolation
method.
"linear" : linear interpolation
"spline" : cubic spline interpolation
"sinc" : standard sin(x)/x interpolation (default)
"sinx" : FFT based periodic sin(x)/sin interpolation
"none" : nearest neighbor (no interpolation)
Returns: A real time series with a PSD identical to the input.
Example:
W1: 1..4
W2: psd(w1)
W3: invpsd(w2)
W4: psd(w3)
W1 contains a 4 point linear ramp. The PSD of the series
is calculated in W2. A time series with the same PSD as
W1 is constructed in W3. The PSD of the constructed time
series is dislayed in W4. Although the values of W1 and W3
differ, the PSD's of the two series are identical.
W2 == {25, 4, 1}
W4 == {25, 4, 1}
Since random phase is used to construct W3, the values will
change if the formula is re-evaluated even if the input does
not change. However, W4 will remain identical to W2 for the
same input.
Example:
W1: 1..4
W2: psd(w1)
W3: invpsd(w2, 0)
W4: psd(w3)
Same as above except the inverse PSD in W3 is calculated with
zero phase and the time series is displayed as an autocorrelation
sequence. The series in W2 and W4 are identical.
W3 == {1.5858, 2, 4.4142, 2}
Since zero phase was used, the values in W3 will remain constant
upon re-evaluation for the same input series. Also notice that
sum(w1) == sum(w3)
Example:
W1: integ(gnorm(1000, 1/1000))
W2: psd(w1)
W3: invpsd(w2, length(w1))
W4: psd(w3);overp(w2);loglog
W1 contains 1000 samples of integrated random noise.
W2 contains the PSD of W1.
W3 returns a time series with the same PSD and length as
the original time series in W1.
W4 compares the PSD of the derived time series with the PSD
of the original time series.
Remarks:
INVPSD finds a time series for the given PSD by reconstructing
the original FFT of the PSD and computing the inverse FFT. Since
the phase information is not available with the PSD, an infinite
number of time series exist with the same PSD.
If N differs from the length of the input time series, the PSD
is interpolated in the frequency domain to produce a time series
of length N.
If rpflag is 1 (the default), a random, odd symmetric phase is
used to reconstruct the time series.
If rpflag is 0, a phase of all zeros is used to reconstruct the
time series and the result is plotted as an autocorrelation
series.
INVPSD always returns a real time series.
See Also:
Fft
Ifft
PSD
#endif
/* calculate an inverse PSD */
invpsd(psd, N, rpflag, method)
{
local iseven, len1, L_odd, L_even, nlen;
local f, df, Fs, s, loff, theta, u, xmin, xmax;
if (argc < 4)
{
if (argc < 3)
{
if (argc < 2)
{
if (argc < 1) error(sprintf("%s - input series required", __FUNC__));
/* time series length */
N = -1;
}
/* random phase */
rpflag = 1;
}
/* interpolation method */
method = "sinc";
}
/* original psd length */
len1 = length(psd);
L_odd = 2 * len1 - 1;
L_even = 2 * len1 - 2;
if (N < 0)
{
/* time series length - default to even */
N = (N < -1) ? bestpow2(L_even) : L_even;
}
if (not(N == L_odd || N == L_even))
{
(xmin, xmax) = xminmax(psd);
if (N % 2 == 0)
{
nlen = N / 2 + 1;
}
else
{
nlen = (N - 1) / 2;
}
f = linspace(0, xmax, nlen);
/* interpolate */
psd = resample(psd, 1/f[2], method);
/* force positive */
psd = abs(psd);
/* delta f */
df = f[2];
setdeltax(psd, df);
}
else
{
df = deltax(psd);
}
iseven = (N % 2 == 0);
/* time series deltax */
Fs = df * N;
/* undo PSD scaling */
s = psd * N * Fs / 2;
/* new psd length */
len1 = length(psd);
/* restore DC and nyquist scaling */
s[1] *= 2.0;
if (iseven)
{
s[len1] *= 2.0;
loff = 2;
}
else
{
loff = 1;
}
/* mag from mag^2 */
s = sqrt(s);
/* form double sided PSD */
s = concat(s, rev(extract(s, 2, len1 - loff)));
if (rpflag)
{
/* add odd symmteric random phase */
theta = grand(int(length(s) / 2), deltax(s), -pi, pi);
if (iseven)
{
theta = {theta, pi, -rev(extract(theta, 2, -1))};
theta[1] = 0;
}
else
{
theta = {0, theta, -rev(theta)};
}
s = s * exp(i * theta);
}
/* inverse transform */
s = real(ifft(s));
/* zero phase - use autocorrelation form */
if (not(rpflag))
{
s = fftshift(s);
}
/* units */
u = extract(psd, 1, 1);
u = sqrt(integ(u));
setvunits(s, getvunits(u));
/* rate */
setrate(s, Fs);
return(s);
}