> > 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?
Yes, this is the Fourier transform of P^|k| (following Ethan D's notation). To derive it you can use the fact that sum(P^k*exp(jwk)) = sum((P*exp(jw))^k) = 1/(1-P*exp(jw)), where the sums are over k >= 0. Then you have to add the same thing for negative frequencies, then subtract a constant because you've accidentally included the k=0 term twice, then turn the algebra crank a couple times. Or get the answer from a table or piece of software. :) 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. 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. 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. Hope that's somewhat helpful! -Ethan On Thu, Nov 5, 2015 at 11:00 AM, Ross Bencina <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 > 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
