>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.
I think you mean to say "infinite energy" and then "energy per unit frequency per unit time," no? E On Thu, Nov 5, 2015 at 8:21 AM, Ethan Fenn <et...@polyspectral.com> wrote: > 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 >
_______________________________________________ dupswapdrop: music-dsp mailing list music-dsp@music.columbia.edu https://lists.columbia.edu/mailman/listinfo/music-dsp
