> > 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. >
Sounds reasonable to me, but I'm not sure I've got it right so who knows! I think I was slightly off when I said that the units of psd are power per unit frequency -- since the whole signal has infinite power, the units really need to be power per unit frequency per unit time, which (confusingly) is the same thing as power. This could be another reason why some special scaling is needed as compared to a finite-length FFT. I'm not sure whether the FFT values should be fringing above the psd line > or not: The psd line is the expected value, so some FFT values should be above it and some below. You could try averaging the squared spectra from a bunch of separate FFT trials and see if that makes things coverge toward the line. -Ethan On Thu, Nov 5, 2015 at 3:48 PM, Ross Bencina <rossb-li...@audiomulch.com> wrote: > 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 >
_______________________________________________ dupswapdrop: music-dsp mailing list music-dsp@music.columbia.edu https://lists.columbia.edu/mailman/listinfo/music-dsp
