Making my own software ASRC (Asynchronous Sample Rate Converter)

[Yansi]

Hello,

Preface:

I am trying to understand, how ASRC (Asynchronous Sample Rate Converter) works. Intention is, to be able to design and write a piece of code, that could solve the very common and typical problem of having a signal crossing two different clock domains with asynchronous clocks. Typically for example, receiving audio from an asynchronously clocked media (receiving it via ETHernet, SPDIF, etc...) and trying to feed them to a DAC, that is getting a local clock, that is of course not synchronous to the source clock.

The usual situation hence looks like I have a receive circular buffer (with DMA) and an output circular buffer (also possibly with DMA). Even though the source and sink sides have a nominal sample rate of say 48 kHz, but the clocks can (and in real world definitely will) drift and have some jitter.  So one of the buffers can potentially overflow (or underflow), depending of which of the clock is faster. I will then be dropping or duplicating samples, which is a big NO as that will result in very audible clicks and signal distortion.

The solution to this according to my research is an ASRC. Even though I have hoarded quite a bunch of different articles, books and other PDFs, I found it difficult to find information about specific implementations (to be able to see examples or more detailed descriptions, other then plain theory.

What I've learnt:

I understand the basic concept of SSRC, when the clocks are synchronous. For example to interpolate by L = 4, I insert three zero samples in between the input samples and than use a band limiting filter to get rid of the in-band image signal. For the band limiting, a FIR filter is typically used. Due to 3/4 of samples being zero, more efficient FIR filter called polyphase is used. I understand the concept of the polyphase FIR filter and how it is possible to reduce the workload by splitting the filter into a L of shorter filters.

What I do not understand fully: .. is how to utilize the polyphase FIR in an ASRC to produce the dynamically variable output length.

What the implementation needs, is a function, that gets fed a constant number of input samples and outputs a variable number of samples (as is required) each time it is called. (Or the other way round: feedin in variable amount of samples and outputting a constant amount each call).
From what I have found, ASRC are being implemented almost the same as SSRC [1], apart that the output polyphase FIR has as many as 2²⁰ phases [2]. But I am not sure, how having a very long output FIR (many phases) allows me to on occasion generate a couple more or less samples. Could anyone please enlighten me how that works? Even though for example the [2] gives some explanation how the coefficient pointer hops in the large table, the explanation seems a bit weird. (See section 3.2 page 18 in [2]).

The only explanation I could come up myself is to have a rather large table of different coefficient sets for the output polyphase FIR (for slightly different ratios) and to switch in between the sets as necessary so, that the amount of output sampes on average will give the required sample rate ratio.

In the end, I will attempt to write down some code and test the ideas above.

[1] lib_src user guide by XMOS, https://www.xmos.ai/download/lib_src-[userguide](1.0.0rc2).pdf
[2] Sample rate conversion in DSPs, Marian Forster, Graz TU Austria, https://www2.spsc.tugraz.at/www-archive/downloads/1_SRC_Marian_Forster_1031275.pdf

[Someone]

What the implementation needs, is a function, that gets fed a constant number of input samples and outputs a variable number of samples (as is required) each time it is called.

At which point you have thrown away the general concept of FPGA processing: streaming processing with minimal memory footprint.
Sorry thought this was in the FPGA forum!

With control of the hardware clocks you can keep a local buffer and change the playback sample rate (slowly) to keep it synchronous in the long term. Slaved to the source material as a PLL.
or
Add ASRC so that the playback sample rate is locked to a local reference, now the same problem as before, you received a chunk of samples from the source but what sample rate are they arriving at to perform the ASRC?

Still need something to extract the actual sample rate of the incoming stream.

[whollender]

ASRCs are tricky.

Let me try to explain it the way I visualize it.

Imagine a poly-phase implementation for a 7/6ths rate change.  Every 6th input sample (and 7th output sample) will be at the same effective sample time, t_0.  Each output sample takes 6/7th of a the input sample time, so the next output sample will precede the next input sample by 1/7th of the input sample period.  Then it will precede by 2/7th, etc, until it lands on the same sample time again.

In this case, each output sample can be at 1 of 7 phases related to the input sample clock, so you need 7 filters that each have slightly different phase to account for the different time delays.

If you then massively scale this up, you can split the input sample period into a large number of 'phases' so that there is a filter for any instantaneous output sample time.  In the reference you linked, they mention 2^20 effective phases.

In a constant rate change system, the way that the output sample phase changes is fixed by whatever the fractional rate change is.  In the above case, the 7/6ths rate change means that each of the 7 phases is used sequentially, but that doesn't need to be the case.  For a 7/3rds rate change, it will start at phase 0, then phase 3/7, then 6/7, then wrap around to 2/7 (9/7), then 5/7, etc.

For an ASRC, there is an additional block that tracks the input and output sample rate, and determines how to move through the 2^20 'phases' that are available.  For instance, if you only had 128 phases, the output phase might tend to jump by 120, so it would be 120, then 112 (120*2 mod 128), then 104, but then it may move by only 119 for a few samples, so it would move to 95 instead of 96, etc.

Because of the scale involved, it's not possible to use normal poly-phase filter techniques.  There are just too many effective phases.  As the second reference describes, with a 64 tap filter, the coefficient ROM for 2^20 phases would require storage for 2^26 coefficients, so they do some interpolation tricks on the fly to generate the exact coefficients from the ROM.

I hope that helps.  I may need to draw some plots to help illustrate the concepts.

In terms of an actual implementation in SW, I don't think you'll want to try to make a function that takes different amounts of input and output data at a time.  I think you'll want to add data to a buffer at a constant rate (the input), then use that input buffer to calculate the filter output at whatever output rate you need.

[Yansi]

Firstly, thank you for responding!

Okay, so the set of FIR filters are in fact fractional delays?  Like in your example, I could have a table of 128 phases. Then what I do, is according to the current ratio of f_out/f_in will choose a phase increment and step through the phases almost exactly as a DDS does. That may start making sense now. If I understood that correctly how that part works.

So, if I am using my brain right:
Lets define the fractional ratio q = f_out/f_in. The phase step then seems to simply be 1/q.

For example:
In case q = 1, the phase increment will be 360° (or 0°). 
In case the ratio will be for example 1.3333  (4/3), i will need to make a 2Pi/q  = 270° step each time.
In case the ratio will be (taken from your example) 7/6, the step phase increment size is 2Pi/q = about 308.6°

Based on the ratio, I calculate the phase increment and pick correct set of coefficients for the FIR. Am I right?
The larger the coefficient table, the smoother the phase stepping. Hence why a specialized DSP chip has the SRC optimized for 2²⁰ phases. (crazy!)

If I have understood the above, then I am missing one more block in the lego: Something to determine the exact ratio q.  But that could be I think determined many ways: Like from the buffer fill rates averaged over time,  direct frequency measurement, etc. I can kind of come up with some ideas.

But back to the polyphase FIR: he big question is: How many phases do I really need to achieve a specific task and how to generate the coefficients?

//EDIT: Like is this just a set of individual fractional delay FIRs [3]? Or is it probably something more complex than that?

[3] Fractional Delay FIR Filters, Neil Robertson, DSPrelated.com, https://www.dsprelated.com/showarticle/1327.php

[Bassman59]

Still need something to extract the actual sample rate of the incoming stream.

Sample rate is easily determined from an AES3 (S/PDIF) stream. Oversample the stream with some high-ish frequency local clock -- say 102.304 MHz, which is 4x the standard 24.576 MHz oscillator and easily generated by an FPGA's PLL or DLL -- and then count the time between edges. The Manchester encoding ensures a transition every bit time. (The high oversampling rate allows this to work for sample rates in the 44.1 kHz family.) Average the counts for the time between transitions and you'll end up with a reasonable estimate of the sample rate.

The larger problem is that there's an implicit assumption that the sample rates are one of the canonical rates: 8 kHz, 16 kHz, 32 kHz and the multiples of 44.1 kHz and 48 kHz. There's no reason one couldn't use a 50 kHz sample rate, but in practice only the "standards" are used. so this simplifies the problem of ASRC to some extent.

[Bassman59]

I've also been interested in implementing ASRC in an FPGA, mainly to avoid using one of the TI SRC4392 chips in a design. Those chips work well, but they're expensive, and if there's already an FPGA in the design, why not use it. The Cirrus parts are half the price but don't have the AES3 receiver built in, but that part too is easily implemented in an FPGA.

The and the Analog Devices [url=https://www.analog.com/en/products/ad1896.html]AD1896 data sheets are good places to start looking at the theory.

Matlab has a feature that generates VHDL an ASRC block. I've looked at it but not implemented it.

[whollender]

Firstly, thank you for responding!

Okay, so the set of FIR filters are in fact fractional delays?  Like in your example, I could have a table of 128 phases. Then what I do, is according to the current ratio of f_out/f_in will choose a phase increment and step through the phases almost exactly as a DDS does. That may start making sense now. If I understood that correctly how that part works.

Yes, conceptually, there are N FIR filters, each with a fractional delay of 2*pi*n/N where n = 0:N-1, and you step through the fractional delays the same way a DDS steps through phase.

So, if I am using my brain right:
Lets define the fractional ratio q = f_out/f_in. The phase step then seems to simply be 1/q.

For example:
In case q = 1, the phase increment will be 360° (or 0°). 
In case the ratio will be for example 1.3333  (4/3), i will need to make a 2Pi/q  = 270° step each time.
In case the ratio will be (taken from your example) 7/6, the step phase increment size is 2Pi/q = about 308.6°

Based on the ratio, I calculate the phase increment and pick correct set of coefficients for the FIR. Am I right?
The larger the coefficient table, the smoother the phase stepping. Hence why a specialized DSP chip has the SRC optimized for 2²⁰ phases. (crazy!)

If I have understood the above, then I am missing one more block in the lego: Something to determine the exact ratio q.  But that could be I think determined many ways: Like from the buffer fill rates averaged over time,  direct frequency measurement, etc. I can kind of come up with some ideas.

Yes, a large part of the design is estimating the actual rate change.  I would look into digital phase-frequency detectors that are used in PLLs to estimate the rate change.

But back to the polyphase FIR: he big question is: How many phases do I really need to achieve a specific task and how to generate the coefficients?

//EDIT: Like is this just a set of individual fractional delay FIRs [3]? Or is it probably something more complex than that?

[3] Fractional Delay FIR Filters, Neil Robertson, DSPrelated.com, https://www.dsprelated.com/showarticle/1327.php

I've never actually needed to implement an ASRC, so you'll probably need to do some experimentation here.  I would probably start with a windowed sinc filter with different fractional delays, ie, using the assumption that the input signal is band-limited to do the interpolation.

With a SW implementation, you could even start by generating the coefficients on the fly for whatever fractional delay you need so that you don't need to actually store a huge coefficient array.

I think most implementations only store some fraction of total number of coefficients and use clever interpolation schemes to generate the coefficients needed for a given phase on the fly.

One thing I forgot in my original post is that this kind of DDS-like fractional delay interpolator does not band-limit the input signal.  If you need to decimate and prevent aliasing, you'll need to filter the input before doing running it through the interpolation filter.  You could also explore the math a bit and try to combine the band-limiting filter and interpolation filter by comparing it to the structure of a poly-phase decimating filter.

The larger problem is that there's an implicit assumption that the sample rates are one of the canonical rates: 8 kHz, 16 kHz, 32 kHz and the multiples of 44.1 kHz and 48 kHz. There's no reason one couldn't use a 50 kHz sample rate, but in practice only the "standards" are used. so this simplifies the problem of ASRC to some extent.

Yeah, it can simplify the design significantly if you know the nominal rates of the input and output, or that the input signal is band-limited at a certain fraction of the input rate.

[Yansi]

Okay, so I will try generating the coefficient sets as the windowed sinc filters, which is exactly what is presented in [3].

Generating the coefficients on the fly seems to be a bit pain, due to the sine function to be calculated per each tap. Now I probably understand, why almost every ASRC implementation I came across was so memory hungry (I mean... couple kB is not a lot, but on a small MCU, it indeed counts!). A lot of coefficient sets and intermediate results to hold!

I've never actually needed to implement an ASRC, so you'll probably need to do some experimentation here.  I would probably start with a windowed sinc filter with different fractional delays, ie, using the assumption that the input signal is band-limited to do the interpolation.

If we assume just cleaning clock jitter/drift of two domains that are supposed to run both at the same nominal rate, is it safe then to just use the polyphase without additional band limiting filter?

From what I understand, the implementation of ASRC in lib-src from XMOS [1] uses dynamically switched filters in front of the "adaptive polyphase".  But that is a wide range ASRC, where that is a must, as the f_out/f_in ratio may be something wild.


I now need to come up with some HW I could put the theory under test and refreshen my matlab (I mean Octave) skills.

[whollender]

Okay, so I will try generating the coefficient sets as the windowed sinc filters, which is exactly what is presented in [3].

Generating the coefficients on the fly seems to be a bit pain, due to the sine function to be calculated per each tap. Now I probably understand, why almost every ASRC implementation I came across was so memory hungry (I mean... couple kB is not a lot, but on a small MCU, it indeed counts!). A lot of coefficient sets and intermediate results to hold!

My suggestion for calculating the coefficients on the fly was mainly to make it easy to experiment with different options while you're still in the software realm.  That should make it easier to iterate different phase steps, different filter lengths, etc.  Also, keep in mind that most HW implementations find ways to reduce the memory footprint by interpolating between filters or coefficients, so you may want to explore that as well in SW to see the effects.

If we assume just cleaning clock jitter/drift of two domains that are supposed to run both at the same nominal rate, is it safe then to just use the polyphase without additional band limiting filter?

Yes, in this case there isn't any out of band signal you need to reject for your new sample rate, so you can just use the polyphase interpolator kernels as is.

From what I understand, the implementation of ASRC in lib-src from XMOS [1] uses dynamically switched filters in front of the "adaptive polyphase".  But that is a wide range ASRC, where that is a must, as the f_out/f_in ratio may be something wild.

I now need to come up with some HW I could put the theory under test and refreshen my matlab (I mean Octave) skills.

Good luck!

[Yansi]

You can start from implementing a linear filter (you can calculate any intermediate points from two points, sample rate and a desired delta T), then filter out the interpolation noise using a heavy FIR filter.

Well, yes, probably a simple good start. What do you mean by a heavy FIR? All I want to do at first is to just make something to be able to compensate drift in between two domains with same nominal wordclock.  I am not sure how one would implement that kind of filter after the lin. interpolator, where the in/out rate is like what, within 500ppm.

If you are decimating, you need to pre-filter the input too.

Yes, of course. But if the ration of interest f_out/f_in is like 0.999?

Polyphase is needed if you need a large sample rate ratio.

Well, if I understand it correctly, what we discussed above is not a polyphase, just an array of different fractional phase delays.

Maybe I should really start with just the linear interpolator. That of course is very easy to understand. But regarding the heavy FIR afterwards, can you elaborate more please, what you meant by that?

[Someone]

Still need something to extract the actual sample rate of the incoming stream.

Sample rate is easily determined from an AES3 (S/PDIF) stream.

The OP specifically mentions an asynchronously arriving source.

You can start from implementing a linear filter (you can calculate any intermediate points from two points, sample rate and a desired delta T), then filter out the interpolation noise using a heavy FIR filter.

Well, yes, probably a simple good start. What do you mean by a heavy FIR?

There are many possible ways to implement the filtering, I think the suggestion is to start with a "toy" example that gets that part of the system working so you can build up all the other parts. Then you can explore what noise/artefacts are produced by different filtering choices.

Higher order FIR (or IIR) filters will better approximate the ideal response but there is no point stressing too much over those details if you don't have anything working.

[Someone]

If we assume just cleaning clock jitter/drift of two domains that are supposed to run both at the same nominal rate, is it safe then to just use the polyphase without additional band limiting filter?

Yes, in this case there isn't any out of band signal you need to reject for your new sample rate, so you can just use the polyphase interpolator kernels as is.

There will be some out of band energy for Fsource > Fsink, but it may be small enough to ignore for the OPs purposes.

[Yansi]

Thanks for all the responses guys! I agree I should start with something very simple, like just the simple linear interpolator.

I need to prepare some hardware to be able to test this easily. I can do some tests in Matlab/Octave, but I am more interested in a solution that can run on real HW. Hold me a thumb, I will try to cobble something fast.

By the way, I thought for a while how to measure the exact in/out wordclock ratio: Not sure if it would work, but how about a FIFO with a PI controller to hold the FIFO half full? The controller will output the ratio to the ASRC.

But not sure if that could work in reality, as there are many unknowns, like for example the time constant of the controller, size of FIFO vs. the amount of of short term jitter and long term drift it can compensate. Everything seems to be some kind of trade-off with these ASRCs.

Guessing the time constant of the controller must be long enough to not distort the audio frequency-wise, like for example second or two, maybe more.

[Yansi]

Would be good also then try to calculate the THD+N for the specific implementation, if thats possible. Like to force matlab/octave to work out the filter and all processing exactly as the MCU will do it, with the limited precision (including the quantization noise). Then analyze the resulting signal for different scenarios.

Otherwise impossible to guess what filtering is enough and what is not.

[gf]

The authors of this paper also give some THD figures for their method which is based on M-fold oversampling, and linear/quadradic/cubic interplation.

[Yansi]

Ah, that paper looks awwsome! Quickly have gone through, it looked full of very usefull stuff. Will read it thoroughly once the HW is done and ready to accept SW.

Hold my beer, tomorrow I have whole day free, so I hope I could come up with some HW.

[gf]

If you are looking for a software solution you could also take a look at various linuxaudio stuff (alsa, jack, pulseaudio, etc.), for instance zita-ajbridge (in conjunction with zita resampler), or Pulseaudio's module-loopback.

[Yansi]

Thank for sharing the links. Will be likely very valuable to study. So far i have just only skimmed through all three. Will get to them back later. I am a bit behind the schedule with making some test HW. I seem to be getting some corona-lockdown-demotivation-syndrome.