I expected this patch to display a constant 0 or maybe a constant 1, but it displays cyclical noise just above 0. Why?
samphold~ noise.pd
-
[samphold~] noise, [phasor~] noise, round-off error, or ?
-
@jameslo because if the period isn't a whole # of samples then the "aliasing" (or time quantization) makes the first sample after the phasor jumps down different every time.
-
I would guess the issue isn't the phasor period.
[metro] can bang only on block boundaries.
If you're using the standard control block size 64 samples, and the typical sample rate 44.1 kHz, then the block does not evenly divide one second. So you're getting snapshots at a period that doesn't exactly equal half a second, and the snapshots would drift relative to the phasor's cycle. (44100 has only 2 prime factors = 2, so the largest power-of-two block size that evenly divides 44100 is 4. Actually I was a bit shocked a few weeks ago to notice that
44100 = 2*2*3*3*5*5*7*7
)But 48000 / 64 = 750 exactly, so that problem might disappear at 48 kHz.
hjh
-
@seb-harmonik.ar But in this case the period is a whole # of samples, 44100, right?
-
@ddw_music I think I sidestepped that block alignment issue by snapshotting the output of [samphold~] at twice the rate it is changing. [samphold~] can't drift relative to the phasor's cycle because it's being triggered by the phasor itself.
-
@jameslo OK, I skipped over that. Must be some roundoff error, then (particularly if phasor~ is adding a phase increment for every sample). I don't have a solid explanation.
hjh
-
@ddw_music If my patch displayed a constant small number, I'd think "Oh, that's right, [phasor~] has to be sub-sample accurate to accommodate non-integral periods, PD could have started it between samples" but that doesn't explain the changing values. At higher frequencies with integral periods, the range of change increases.
-
@whale-av said:
@jameslo It's not a rounding error........ it is drifting for some reason.
https://github.com/pure-data/pure-data/blob/d07e6efbaef57442d11ca6a83810187256db6292/src/d_osc.c#L95
... seems to show that phasor~ is integrating floating point values.
These floating point values necessarily have some inaccuracy. So their integration will have some inaccuracy, which will accumulate.
If the result of adding 1/44100 to itself 44100 times is not exactly 1.0, over time this would manifest as drift.
That is, rounding error can be drift. ...?
hjh
-
@jameslo of course, whoops. in that case my guess would agree w/hjh, probably some case of the phase increment/state not being able to be represented exactly internally
-
@ddw_music Maybe, but it is the same calculation every time so why would the errors accumulate?
No problem with [line~] or [vline~] into [samphold~]
[phasor~] phase not updated except at block boundary? I tried a single sample block and there was no difference........
Changing the Pd samplerate does compress the pattern. At 200K sample rate the deviation is considerably reduced.
FFT...... I agree..... still probably phase.
Does anyone have [vphasor~] compiled?
https://forum.pdpatchrepo.info/topic/10192/vphasor-and-vphasor2-subsample-accurate-phasors
David. -
@whale-av I think @ddw_music is saying that each point in the phasor except for the jump downward is generated by adding an incremental amount to a running subtotal. That amount is a function of the frequency, and if its calculation contains round-off error, then that error accumulates in the subtotal. It's consistent with that nice graph you showed, assuming that the jump downward is produced by something that doesn't zero out the accumulated error.
I think this patch shows that [phasor~ 1] drops a sample every 5 seconds or so, i.e. it runs a little fast:
phasor~ timing.pd
I suppose I could/should get off my lazy rear, learn something about github, download the source, and search through it, but I'd rather continue slogging through katya's blog. -
@jameslo the phasor~ source was a bit confusing to me at first, it works by setting the upper 32 bits of a double (including the exponent) such that bit 32 of the double represents the "1s" place value, so all of the lower 32 bits represent the fractional part. Then after the correct phase accumulation is added, the magnitude that was added to make all lower 32 bits the fractional part is subtracted from the double in order to output just the fractional part. (and then the upper 32 bits are again set to that value, destroying the non-fractional part)
All of this is done to avoid cpu branching so it can pipeline instructions better (so instead of having to check if the phase > 1, each sample you just set bits to a certain value) -
@jameslo just investigating further: the closest that a 32 float value can come to representing 1/44100 seems to be 0.00002267573654535226523876190185546875. However, since bit value 32 must be 1 and because the exponent of this number is 2^-16, that means that once this number is added to 1572864 as a double the fractional part can only represent
0.00002267584204673767089843750.0000226758420467376708984375×44100 = 0.00000463426113128662109375
as the fractional part of the first "wrapped-around" sample (I might be missing something, but something like that is probably happening)
edit: actually after testing this matches exactly with the behavior of your first patch
(1st wrap-around is that # after sending "1" to the phase of phasor~ and the toggle at the same time) -
@seb-harmonik.ar (full disclosure: I'm barely hanging on...) OK, is it consistent with my last patch, the one where the phasor falls short one sample every 5 seconds? I didn't say it, but that seemed like way more than accumulated single-precision floating point round-off error to me, which is why I'm insecure about my conclusion.
Just to clarify, are you saying that the accumulator is single or double precision floating point? Maybe even simpler and more to the point, are you confirming @ddw_music's round-off error accumulation theory?
-
@jameslo yes, @ddw_music was correct. It's a combination of double and single precision (and I edited my last comment, this prediction is consistent w/ the 1st patch at least). The value "conv" in the source code is stored in a t_float, which is generally a 32-bit float these days (though maybe that will change soon..). This is set to 1/samplerate. every sample conv is multiplied by the input frequency (which is also a 32-bit float) and then added as a double to the current phase, which is a double with value 1572864 + actual phase, (1572864 is 3^19, a float value that makes bit 32 of the entire 64-bit double value the 1s place, leaving the remaining 32 bits as the fractional part). Every sample the top 32 bits of the phase are set to the top 32 bits of 1572864, and when the phase is output, 1572864 is subtracted from it.
tldr: the phase accumulator is basically 32-bit fixed-point
-
@jameslo Just for information.
If the sample rate is a power of 2...... 32768Hz for example...... there is no drift.
I wonder whether your RMS an FFT results would converge at such a sample rate?
David. -
@whale-av that makes sense bc powers of 2 are easily represented exactly in binary floating point (so the increment can be represented exactly without roundoff error)
-
@whale-av said:
I wonder whether your RMS an FFT results would converge at such a sample rate?
He he he, they don't. I'm glad because if they did my head would explode.
@seb-harmonik.ar said:
the closest that a 32 float value can come to representing 1/44100 seems to be 0.00002267573654535226523876190185546875. However, since bit value 32 must be 1 and because the exponent of this number is 2^-16, that means that once this number is added to 1572864 as a double the fractional part can only represent
0.0000226758420467376708984375Wait, that is consistent with my second patch--the increment is slightly large, and so the phasor reaches its peak a bit early, i.e. its running a bit fast. In numbers, the increment is ~ 0.106e-9 too big. That times 44100 times 5 is greater than 23e-6, which is greater than 1/44100 as a 32 float, i.e. 1 sample.
