Post B6dAfI3hOLzRnDRnXc by azonenberg@ioc.exchange
(DIR) More posts by azonenberg@ioc.exchange
(DIR) Post #B6cvEQoZEN0dol6X56 by azonenberg@ioc.exchange
0 likes, 0 repeats
Fun engineering problem I am currently facing: Generate a very long sine wave in ngscopeclient, accurately and in parallel, using only single-precision floating point and 32-bit integer math.The current implementation, used in both the "sine" and "downconvert" filter, looks basically like (pseudocode with some stuff like scaling and offsets removed for simplicity)double rate;uint depth;for GPU thread i in 0...depth-1 double tmp = (i * rate) mod 2*pi out[i] = sin( (float) tmp)The problem is the use of the float64 intermediate, which breaks on platforms (most notably apple silicon) that do not have native fp64 support.Naively having rate and tmp be float32 causes loss-of-precision errors and waveform distortion with deep memory.Anybody have ideas on some clever numerical method to do this with only float32 temporaries?
(DIR) Post #B6cvcEABv3j3WxZ84G by azonenberg@ioc.exchange
0 likes, 0 repeats
If this were being done sequentially, it would be straightforward to do something likefloat phase = 0for i in 0...depth-1 out[i] = phase kahanSum(phase, rate) if(phase > 2*pi) phase -= 2*pibut I'm not sure how this would translate to a parallel implementation
(DIR) Post #B6cvfpwnta77aJ1bIe by mattieuMattieu@mstdn.social
0 likes, 0 repeats
@azonenberg this is the classic technique and is still very useful. https://ir.cwi.nl/pub/9159/9159D.pdf
(DIR) Post #B6cvtUdbv929Y3YjTc by dlharmon@chaos.social
0 likes, 0 repeats
@azonenberg For N instances, do starting phases of rate*i, increment phase by N*rate in each instance?For long running sine generation, I always use integers for phase, scale such that integer overflow does the mod 2*pi for me.
(DIR) Post #B6cwJs5TMJoDnwAZOq by azonenberg@ioc.exchange
0 likes, 0 repeats
Just to illustrate the difference in precision, here's the very end of a 100 million sample sinewave with float32 and float64 temporaries.
(DIR) Post #B6cwcbna8xwEmwjgq8 by jawnsy@mastodon.social
0 likes, 0 repeats
@azonenberg for auld long sine
(DIR) Post #B6cwj9DSoMHkAUxZuS by pkhuong@discuss.systems
0 likes, 0 repeats
@azonenberg Rate * large integer is already a problem, you want to commute with the range reduction. Random idea: precompute rate * (2**k) mod (2 pi) on the host, at high precision, then sum up to (log n) LUT entries before another range reduction.
(DIR) Post #B6cwq2SSKuzfEeG53o by azonenberg@ioc.exchange
0 likes, 0 repeats
@dalias Max record length is 4GB with the current programming model so 2^30 float32 samples. I have so far not seen a problematic loss of precision with float64 temporaries using this algorithm, I just didn't realize some of our target GPUs did not have float64 support until recently.
(DIR) Post #B6cx52ucvg7KzKai5Q by azonenberg@ioc.exchange
0 likes, 0 repeats
@dalias A small static frequency error is better than the output not being sinusoidal. But the frequency changing over the duration of the waveform isn't great.
(DIR) Post #B6cxZIR3mN9XotBaSG by datenwolf@chaos.social
0 likes, 0 repeats
@azonenberg Two years ago I had to solve a similar problem. Implementation details are subject of my PhD thesis. Cliffs notes version: instead of a rate think in terms of fractional interval lengths. Then map onto a skewed 2D grid, with truncated stride, and compensation for kerfed samples in calculation of sample address. 2D grid maps to in-period samples / periods.
(DIR) Post #B6d1pLfduzBCIsutY8 by azonenberg@ioc.exchange
0 likes, 0 repeats
So I ended up doing a somewhat simple solution: divide the index into blocks of N (currently 32768 but I can play with this) and then calculate the fraction of a cycle, from 0 to 1, within each block.Then I can calculate the final phase as a*block_index + b*position_within_block where all of the values are relatively small and rounding error isn't significant. Then multiply by 2*pi and sine it.This seems to be sufficient for now given the 2^30 maximum input record length. It's not worse than what I had before with float64, and works entirely with float32 temporaries.
(DIR) Post #B6d2AVirjxKDYoVziy by ChuckMcManis@chaos.social
0 likes, 0 repeats
@azonenberg I was trying to figure out what you were optimizing for, phase precision? frequency precision? minimal quantization harmonics? Something else?
(DIR) Post #B6d2YYVz0hwjAcxG9A by ponygol@chaos.social
0 likes, 0 repeats
@azonenberg Hard limit on float32 or are 32-Bit integer types also fine? Seems like calculating tmp in 64-Bit FXP using two 32-Bit uint32 should be possible at negligible performance loss.
(DIR) Post #B6d2deXg1p3xJOB4eu by azonenberg@ioc.exchange
0 likes, 0 repeats
@ponygol f32 and i32 are universally available f64 and u64 are optional extensions
(DIR) Post #B6d2lqCeYc5fnIFSPw by azonenberg@ioc.exchange
0 likes, 0 repeats
@ChuckMcManis smooth sine with minimal distortion and discontinuities is the main goal but not having significant frequency error is also important
(DIR) Post #B6d3HszvhZYNjEZqwy by ponygol@chaos.social
0 likes, 0 repeats
@azonenberg If you get the rate in Hz instead of radial frequency, you can replace the modulus by just taking the fractional part.If you convert rate (which is most likely a lot less than 1) into Q0.31, you can then just multiple by i and ignore all overflow (this is the modulus). Convert the result back to float, multiply by 2*pi and calculate the sine?
(DIR) Post #B6d3ovX3Igvz91pDua by ponygol@chaos.social
0 likes, 0 repeats
@azonenberg (there might be some care needed to see how negative numbers are handled during overflow, but I currently am not in a place to work out the details)
(DIR) Post #B6d3ovi2dpizh6y02a by azonenberg@ioc.exchange
0 likes, 0 repeats
@ponygol i should be able to do it all in u32 domain with no negative values ever. Will try in a bit busy with other stuff rn
(DIR) Post #B6d4EgN0IMo1eKPFT6 by AMS@infosec.exchange
0 likes, 0 repeats
@azonenberg Usually I do this with sine/cosine taking periods instead of radians, but that's immaterial. You'll want to split rate and i into two half-precision parts (12 mantissa bits each iirc) then tmp= i*rate becomes tmp_high = i_high * rate_high + i_high*rate_low + i_low*rate_hightmp_low = i_high * rate_high - tmp_high + i_high*rate_low + i_low*rate_high + i_low*rate_lowtmp = (tmp_low mod 2pi) + (tmp_high mod 2pi)P.S.: Yes, I realize this is "we have fp64 at home".
(DIR) Post #B6d6IanR9xojQhlKd6 by azonenberg@ioc.exchange
0 likes, 0 repeats
@AMS I'm no stranger to shader based software bighum for int64 lol
(DIR) Post #B6dAAedpNbwWr3jp4K by penguin42@mastodon.org.uk
0 likes, 0 repeats
@azonenberg Can you do something like bresenhams algorithm?
(DIR) Post #B6dAe13H8Zi8l1dqPg by azonenberg@ioc.exchange
0 likes, 0 repeats
@penguin42 bresenham is iterative so you can't jump to the nth sample in parallel random access
(DIR) Post #B6dAfI3hOLzRnDRnXc by azonenberg@ioc.exchange
0 likes, 0 repeats
@penguin42 this would be much easier sequential
(DIR) Post #B6dCngXcRwyPwuZajY by ChuckMcManis@chaos.social
0 likes, 0 repeats
@azonenberg Okay. Is 'Frequency Error': |desired_f - actual_f|In spectrum terms, is 'minimal distortion' equal to the Spurious-Free Dynamic Range (SFDR)? (trying to map your requirements into things I've worked on in my Quadratic NCO)
(DIR) Post #B6dEN1vN2WrcMD6H7g by azonenberg@ioc.exchange
0 likes, 0 repeats
@ponygol yep, this worked out pretty nicely//Get the fractional phase using integer math with implicit mod 2^32 to handle wrappinguint fpfrac = i * fpfreq;const float fix_to_float = 1.0 / 4294967295.0;float frac = float(fpfrac) * fix_to_float;//Now that we have the fractional phase,//computing the sine can be done in the float32 domain without loss of precisionfloat two_pi = 6.28318530717;dout[i] = bias + (scale * sin(frac*two_pi + startphase));
(DIR) Post #B6f25geaXSnvh4su0m by encthenet@flyovercountry.social
0 likes, 0 repeats
@azonenbergI'm not familiar with floating point precision, but I take it that doing(i mod 2*pi) * (rate mod 2*pi)Doesn't solve things?
(DIR) Post #B6fEOUgyQzynY1jbma by jenesuispersonne@piaille.fr
0 likes, 0 repeats
@azonenberg I don't know if it may help but projects like GNURadio or VCVrack (or any synth. modeler) should also face same kind of problems to minimize signal processing errors.
(DIR) Post #B6fGLJOYBFGPY9jAZs by azonenberg@ioc.exchange
0 likes, 0 repeats
@jenesuispersonne I went with the idea @dlharmon suggested, having the phase accumulator be fixed point 0.32 resolution. This lets you do all the intermediate "get phase for each sample" math with uint32s, no loss of precision, and the mod 2pi happening naturally with 32-bit overflow.Then just rescale to float32, multiply by 2pi over 2^32-1 , and you're golden.