Yes, thank you! I guess most of the places I typed the word power I really meant energy... units are hard...
-Ethan On Thu, Nov 5, 2015 at 7:33 PM, Ethan Duni <ethan.d...@gmail.com> wrote: > >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 >
_______________________________________________ dupswapdrop: music-dsp mailing list music-dsp@music.columbia.edu https://lists.columbia.edu/mailman/listinfo/music-dsp
