Thanks Ethan,
I think that I have it working. It would be great is someone could check
the scaling though. I'm not sure whether the FFT values should be
fringing above the psd line or not:
https://www.dropbox.com/s/txc0txhxqr1t274/SH1_2.png?dl=0
I removed the hamming window, which was causing scaling problems. The
FFT output is now scaled so that the the sum of power over all bins
matches the power of the time domain signal:
https://gist.github.com/RossBencina/a15a696adf0232c73a55/bdefe5ab0b5c218a966bd6a04d9d998a708faf99
On 6/11/2015 12:02 AM, Ethan Fenn wrote:
And is psd[w] in exactly the same units as the magnitude squared
spectrum of x[n] (i.e. |ft(x)|^2)?
More or less, with the proviso that you have to be careful whether
you're talking about power per unit frequency which the psd will give
you, and power per frequency bin which is often the correct
interpretation of magnitude squared FFT results -- the latter depending
on the FFT scaling conventions used.
Let's see if I got this right: each bin contains the power for a
frequency interval of 2pi/N radians. If I multiply each bin's power by
N/2pi I should get power values in units of power/radian.
The psd makes no reference to any transform length, since it's based on
the statistical properties of the process. So I think it would be wrong
(or at least inexact) to have a scale related to N applied to it. If you
want the magnitude squared results of an FFT to match the psd, it seems
more correct to scale the FFT and try a few different N's to see what
factor of N will give consistent results.
That makes sense.
As to the exact scale that should be applied... I think there should be
a 1/3 in the expression for psd, because E[x^2]=1/3 where x is uniform
in [-1,1]. Aside from that, there might be a factor of 2pi depending on
whether we're talking about power per linear or angular frequency. And
there could be others I'm not thinking of.... maybe someone else can
shed more light here.
I multiplied the psd by 1/3 and as you can see from the graph it looks
as though the FFT and the psd are more-or-less aligned.
Hope that's somewhat helpful!
Very clear thanks,
Ross.
-Ethan
On Thu, Nov 5, 2015 at 11:00 AM, Ross Bencina
<rossb-li...@audiomulch.com <mailto:rossb-li...@audiomulch.com>> wrote:
Thanks Ethan(s),
I was able to follow your derivation. A few questions:
On 4/11/2015 7:07 PM, Ethan Duni wrote:
It's pretty straightforward to derive the autocorrelation and
psd for
this one. Let me restate it with some convenient notation. Let's say
there are a parameter P in (0,1) and 3 random processes:
r[n] i.i.d. ~U(0,1)
y[n] i.i.d. ~(some distribution with at least first and second
moments
finite)
x[n] = (r[n]<P)?x[n-1]:y[n]
Note that I've switched the probability of holding to P from
(1-P), and
that the signal being sampled-and-held can have an arbitrary (if
well
behaved) distribution. Let's also assume wlog that E[y[n]y[n]] = 1
(Scale the final results by the power of whatever distribution
you prefer).
So for y[n] ~U(-1,1) I should multiply psd[w] by what exactly?
Now, the autocorrelation function is ac[k] = E[x[n]x[n-k]].
Let's work
through the first few values:
k=0:
ac[0] = E[x[n]x[n]] = E[y[n]y[n]] = 1
k=1:
ac[1] = E[x[n]x[n-1]] = P*E[x[n-1]x[n-1]] + (1-P)*E[x[n-1]y[n]] =
P*E[y[n]y[n]] = P
The idea is that P of the time, x[n] = x[n-1] (resulting in the
first
term) and (1-P) of the time, x[n] is a new, uncorrelated sample from
y[n]. So we're left with P times the power (assumed to be 1 above).
k=2:
ac[2] = P*P*E[x[n-2]x[n-2]] = P^2
Again, we decompose the expected value into the case where x[n]
= x[n-2]
- this only happens if both of the previous samples were held
(probability P^2). The rest of the time - if there was at least one
sample event - we have uncorrelated variables and the term drops
out.
So, by induction and symmetry, we conclude:
>
ac[k] = P^abs(k)
>
And so the psd is given by:
psd[w] = (1 - P^2)/(1 - 2Pcos(w) + P^2)
What is the method that you used to go from ac[k] to psd[w]? Robert
mentioned that psd was the Fourier transform of ac. Is this
particular case a standard transform that you knew off the top of
your head?
And is psd[w] in exactly the same units as the magnitude squared
spectrum of x[n] (i.e. |ft(x)|^2)?
Unless I've screwed up somewhere?
A quick simulation suggests that it might be okay:
https://www.dropbox.com/home/Public?preview=SH1_1.png
But I don't seem to have the scale factors correct. The psd has
significantly smaller magnitude than the fft.
Here's the numpy code I used (also pasted below).
https://gist.github.com/RossBencina/a15a696adf0232c73a55
The FFT output is scaled by (2.0/N) prior to computing the magnitude
squared spectrum.
I have also scaled the PSD by (2.0/N). That doesn't seem quite right
to me for two reasons: (1) the scale factor is applied to the linear
FFT, but to the mag squared PSD and (2) I don't have the 1/3 factor
anywhere.
Any thoughts on what I'm doing wrong?
Thanks,
Ross.
P.S. Pasting the numpy code below:
# ---------------8<--------------------------------
# see
https://lists.columbia.edu/pipermail/music-dsp/2015-November/000424.html
# psd derivation due to Ethan Duni
import numpy as np
from numpy.fft import fft, fftfreq
import matplotlib.pyplot as plt
N = 16384*2*2*2*2*2 # FFT size
y = (np.random.random(N) * 2) - 1 # ~U(-1,1)
r = np.random.random(N) # ~U(0,1)
x = np.empty(N) # (r[n]<P)?x[n-1]:y[n]
# hold process with probability of P for holding
P = 0.9
x[0] = 0
for i in range(1,N-1):
if r[i] < P:
x[i] = x[i-1]
else:
x[i] = y[i]
# FFT of x
window = np.hamming(N)
spectrum = fft(x*window,N)
magnitudeSq = ((2.0/N) * np.abs(spectrum)) ** 2
# Calculate PSD
w = fftfreq(N)*np.pi*2 # (w from 0 to 0.5 then -0.5 to 1
psd = ((1.0 - P**2)/(1.0 - 2*P*np.cos(w) + P**2))
psd = psd*(2.0/N) # not sure about this
plt.plot(w, magnitudeSq)
plt.plot(w, psd, color='red')
plt.xscale('log')
plt.yscale('log')
plt.show()
# ---------------8<--------------------------------
_______________________________________________
dupswapdrop: music-dsp mailing list
music-dsp@music.columbia.edu <mailto:music-dsp@music.columbia.edu>
https://lists.columbia.edu/mailman/listinfo/music-dsp
_______________________________________________
dupswapdrop: music-dsp mailing list
music-dsp@music.columbia.edu
https://lists.columbia.edu/mailman/listinfo/music-dsp
_______________________________________________
dupswapdrop: music-dsp mailing list
music-dsp@music.columbia.edu
https://lists.columbia.edu/mailman/listinfo/music-dsp
