Old 12-23-2018, 11:17 PM   #1
clepsydrae
Human being with feelings
 
clepsydrae's Avatar
 
Join Date: Nov 2011
Posts: 2,236
Default upsampling as interpolation: what LPF?

TLDR: What's a smart brickwall filter to use for LPF after upsampling which I can use in JSFX without an advanced degree in math?

-----

I want to upsample an array of data representing a continuous signal to find the "true intersample peak". So I zero-stuff the data into a buffer and low pass filter it. It sounded so simple on paper. :-)

(Such data could be audio, but in this case it happens to be the results of autocorrelation for pitch detection; range is not +/- 1).

I'm finding the filter coefficients in python (scipy.signal), and copy/pasting in to JSFX, so I'm sure I'm losing some precision there, and I'm wondering if that explains my issues.

The data set is so small (on the order of 50-1000 points after upsampling) that the filter can't seem to settle in usefully. some filters settle in well, but only for some frequencies, etc.

I've been trying all manner of IIR and FIR filters, tweaking all their parameters, etc. They work fine for regular audio at regular sample rates, albeit with some ringing artifacts. But on the autocorrelation data set they don't work well: tons of ringing, or too slow to settle, etc. I started with steep elliptic filters near nyquist with little luck. My best effort now is a butterworth at 1/8 nyquist. The data is pretty smooth already, so I can afford to lose HF, but I don't really get why I need to... Brickwall LPF after upsampling happens all the time, and it can't be as destructive as what I'm seeing... What's the magic trick? :-) Is there some bog standard?

I can't get FIR filters to work at all on the data: they generate regularly spikey results at the same interval as the filter delay line, or other weird results.

My not-really-success is happening at like 4X upsampling. I'd like to go much higher, but everything falls apart when I try.

In short: what's the best way to interpolate a short data set to smooth it out?

If my problem is indeed the copy/paste from python, would upsampling/filtering in several stages be smarter?

This is the simple pretty little picture I'm trying to reproduce. What voodoo is this?

clepsydrae is offline   Reply With Quote
Old 12-24-2018, 04:00 AM   #2
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Germany
Posts: 6,312
Default

As upsampling is a well known mathematical / scientific issue, I am sure that there is decent literature about it, and I don't suppose there is any matter of "taste" with it but a single correct mathematical / scientific solution, as there is only one solution for the equation resulting in the original band limited signal from the samples.

If band limiting is done correctly, the A/D <-> D/A process is reversible (and hence unique).

As this is true for any sample rate, you don't have any choice with the upsampling procedure. And Upsampling will not help in any way, unless the plugin or hardware to be fed with the upampled signal is not able to decently work with the lower sample rate.

-Michael
__________________
www.boa-sorte.de
mschnell is online now   Reply With Quote
Old 12-24-2018, 04:35 AM   #3
sai'ke
Human being with feelings
 
sai'ke's Avatar
 
Join Date: Aug 2009
Location: Germoney
Posts: 195
Default

Re the previous answer, there are always design constraints. Ideal interpolation would be with an infinitely long sinc, which is not implementable in software or real-time. You can approximate it somewhat with a filter, or do stuff with zeroing parts of the FFT, but then there's a huge latency involved. I mean, sure, if you know exactly how bandlimited the signal is, then you could design a filter that does near ideal filtering, but I don't think it's ever fully reversible.

As for upsampling, either FIR or IIR should work. I am skeptical whether loss of precision due to c/p-ing filtercoefficients is much of a factor. For the FIR (and maybe also the IIR), I suspect that there is a bug in your filter code.

What I would suggest is building the filter code you plan to use in JSFX in python first and checking it against filtering with python's filter library in python, just to make sure that part of the implementation is correct and don't write things up to numerical precision too quickly. Having a ground truth is extremely helpful, so doing it in python first might save you time rather than cost more

That said, I would always start with the simplest approach, is linear interpolation or quadratic not fine enough? I mean, sure you can probably get more accuracy out of a more ideal filter, but also more "weird things that can happen".

You might be interested in this too: https://forum.cockos.com/showthread.php?t=175038
I think both Tale and SaulT have upsampling libraries with pretty steep filters. You might be able to find more threads of either of them that show some more details. I also implemented upsampling in my Filther. I remember playing with steeper filters before, but ending up settling with a cheaper filter based on graphdist.
__________________
[Soundcloud] | [Tracker Plugin: Thread|Github|Reapack link] | [Machine UI Plugin: Thread|Github|Reapack link] | [Filther: Github]

Last edited by sai'ke; 12-24-2018 at 05:14 AM.
sai'ke is offline   Reply With Quote
Old 12-24-2018, 05:22 PM   #4
clepsydrae
Human being with feelings
 
clepsydrae's Avatar
 
Join Date: Nov 2011
Posts: 2,236
Default

I think I'm designing the filters wrong... and/or the precision issue isn't helping...

Quote:
What I would suggest is building the filter code you plan to use in JSFX in python first and checking it against filtering with python's filter library in python, just to make sure that part of the implementation is correct and don't write things up to numerical precision too quickly.
Wise advice, thanks, hopefully I can get to that.

So, the filtering does work fine in JSFX already, as long as I don't upsample. That is, I can apply these FIR and IIR filters at whatever frequency, whatever number of poles, etc, and they work as expected. At 2X upsampling, the result gets a little fuzzy, and it falls apart at 4X and 8X.

Examples: original in blue, upsampled in red; these results from a 15-order FIR linear phase low pass at 20k with a 2k transition; these are the "best" results I've generated, most are crazy --

1X -- no upsampling, just the filter:

2X upsampling, filter re-designed with appropriately scaled cutoff, i.e. f/2:

4X upsampling:

8X upsampling:


I've tried this with higher/lower order, with IIR, FIR, different frequency ranges, etc.

Quote:
That said, I would always start with the simplest approach, is linear interpolation or quadratic not fine enough? I mean, sure you can probably get more accuracy out of a more ideal filter, but also more "weird things that can happen".
I don't think linear would work for me as I'm seeking the X location of the signal peak, and that would just return the same result, yeah? I've considered simpler methods, but at this point I need to figure this out for my own education and sanity. :-)

Quote:
You might be interested in this too: https://forum.cockos.com/showthread.php?t=175038
D'oh! Can't believe I didn't find that already! Works great. I'm upsampling 64X and it's nice and smooth. I'm going to go forward with it but I'd still love to know where I'm going wrong with my testing.

Question: the image shown in that thread shows the rolloff starting around 20k and ending around 40k. I thought the lowpass in upsampling needed to cut off before nyquist to filter out mirrored spectra. E.g. for 44.1k you have to remove everything 22050Hz and above? If that image doesn't assume 44.1k sample rate, then where is nyquist for that graph?

Question 2: whether my own test or the library you mentioned, I'm manually aligning the filtered samples with the original to correct for the filter delay. Is there some way to know this figure exactly? For 64X, using the library, it's ~566 samples... I'm not seeing an obvious way to derive that from 64 or the filter order (19, I believe). I'll need to know the exact amount for my purposes...

Thanks in advance!
clepsydrae is offline   Reply With Quote
Old 12-25-2018, 04:14 AM   #5
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Germany
Posts: 6,312
Default

Quote:
Originally Posted by sai'ke View Post
Re the previous answer, there are always design constraints. Ideal interpolation would be with an infinitely long sinc, which is not implementable in software or real-time. You can approximate it somewhat with a filter, or do stuff with zeroing
I suppose for offline application, a sinc long enough to rather correctly reproduce the amplitude resolution (e.g. 24 bits) should be the way to go.

Is a sic approximation in fact a normal (low pass) "filter", (FIR / IIR) mathematically ?

Quote:
Originally Posted by sai'ke View Post
That said, I would always start with the simplest approach...
If there is a decent way to judge the result.

Quote:
Originally Posted by sai'ke View Post
is linear interpolation or quadratic not fine enough?
A very usual (non-audio) interpolation is cubic spline which creates good optical results with not too much calculation overhead. But I don't know if this really works well for audio.

-Michael
__________________
www.boa-sorte.de

Last edited by mschnell; 12-25-2018 at 04:26 AM.
mschnell is online now   Reply With Quote
Old 12-25-2018, 04:42 AM   #6
sai'ke
Human being with feelings
 
sai'ke's Avatar
 
Join Date: Aug 2009
Location: Germoney
Posts: 195
Default

@mschnell
Yeah, a sinc is a lowpass filter. An infinite sinc would be the ideal brick-wall lowpass filter (perfect cutting). Usually people window it, which makes it slightly less ideal.

@clepsydrae
Ah yes, sorry, I thought you were looking for zero crossings . If you're looking for peaks, linear interpolation makes zero sense.

I suspect that phase is important to you in this case, so I would go with a FIR filter. A FIR filter has linear phase and the group delay is ( N − 1 ) / 2 with N being the number of taps. They're a little less effective at cutting, but that's the price you pay. You can even do the upsampling and filtering in a single step with a FIR.

Alternatively, if it's not realtime, what one can do is filter with an IIR, then reverse the signal and filter again. This gets rid of any group delay and gives you twice the filter effect (also transients on both ends tho). Or use the FFT. FFT the signal, zero pad it, and then transform back.

I can't tell what's going wrong just from those graphs.
__________________
[Soundcloud] | [Tracker Plugin: Thread|Github|Reapack link] | [Machine UI Plugin: Thread|Github|Reapack link] | [Filther: Github]
sai'ke is offline   Reply With Quote
Old 12-25-2018, 02:46 PM   #7
clepsydrae
Human being with feelings
 
clepsydrae's Avatar
 
Join Date: Nov 2011
Posts: 2,236
Default

Quote:
Originally Posted by sai'ke View Post
A FIR filter has linear phase and the group delay is ( N − 1 ) / 2 with N being the number of taps.
Thanks! I forgot that the upsampling library requires recursive 2X cascades, so I need to consider the delay at each stage.

Quote:
I can't tell what's going wrong just from those graphs.
I don't know how much time you have to devote to my education, but in case you have an inspiration to see where I'm going wrong, here are a couple demos:

UpsampleDemo.txt -- upsamples 8X, applies the filter, draws the original and the result (what I used for the images I posted.)

FIRDemo.txt -- uses the same filter, but the coefficients generated for a non-upsampled sample rate. Judging by white noise through it, it seems to work as expected.

That's not the greatest filter, but again, I've tried 200 combinations of parameters, and maybe I've just missed the right one but I can't seem to get anything to work.

For the record, those coefficients were generated with scipy's signal.firwin:

Code:
from scipy import signal
numtaps=25
cutoff=0.5442176870748299  # 12k at 44.1k sample rate
width=0.1360544217687075  # 3k transition at 44.1k sample rate
upsamp=1

b = signal.firwin(numtaps, cutoff/upsamp, width/upsamp)

for i in range(len(b)): print('filtb[' + str(i) + '] = {:.100f}'.format(b[i]) + ';')
clepsydrae is offline   Reply With Quote
Old 12-25-2018, 09:39 PM   #8
clepsydrae
Human being with feelings
 
clepsydrae's Avatar
 
Join Date: Nov 2011
Posts: 2,236
Default

Also a note for anyone finding this thread: I learned that you may want to normalize your data before upsampling (perhaps pre-setting the filter taps to the starting value of your data set would solve this). Mine had a large range (e.g. min 1400 to max 1600) and the results were fuzzy and had a strong ring at the start -- I normalized to +/- 1 and all is well.

Last edited by clepsydrae; 12-25-2018 at 11:02 PM.
clepsydrae is offline   Reply With Quote
Old 12-26-2018, 01:38 AM   #9
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Germany
Posts: 6,312
Default

Quote:
Originally Posted by sai'ke View Post
@mschnell
Yeah, a sinc is a lowpass filter. An infinite sinc would be the ideal brick-wall lowpass filter (perfect cutting). Usually people window it, which makes it slightly less ideal.
Obviously.
Can a decent "inaccuracy" be calculated from the window size to allow for selecting the best quality / CPU overhead ratio ?
For a test I once did a sinc filter with FFT in JSFX. Here the FFT Window just sets the low limit of the spectrum that is handled.
Are there ways to do a sinc filter without doing a Fourier transform ?

-Michael
__________________
www.boa-sorte.de

Last edited by mschnell; 12-26-2018 at 03:40 AM.
mschnell is online now   Reply With Quote
Old 12-26-2018, 03:36 AM   #10
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Germany
Posts: 6,312
Default

Doing more reasoning...

I am not convinced regarding the sinc (or any) filter solution.

Due to theory the samples define a unique band limited waveform and same is identical to that band limited waveform which had been sampled.

To reproduce the waveform you need to feed a band limiting "thing" with (as short as possible) pulses with the "power" of any sample.

The times between the pulses would be "open" (electrically spoken) or "missing value" (mathematically spoken).

To use any kind of "normal" (electronic or software) filter, you would need to provide values for the time between two samples ("hold", linear interpolation, quadratic interpolation, spline interpolation, ...) this obviously creates a distortion to the signal to be reproduced, even if perfect band-limiting might keep the damage at bay. (See the "adding Zeros" picture in the first post here, which does not state what exactly the "interpolation filter" is supposed to be. But maybe adding zeros can be proven to cause a gain reduction as the only "distortion".)

I suppose there should be an algorithm that can cope with missing values. I understand that something like this is used in "oversampling D/A converters".

-Michael
__________________
www.boa-sorte.de

Last edited by mschnell; 12-26-2018 at 11:36 AM.
mschnell is online now   Reply With Quote
Old 12-27-2018, 01:17 AM   #11
clepsydrae
Human being with feelings
 
clepsydrae's Avatar
 
Join Date: Nov 2011
Posts: 2,236
Default

Son of a gun, I found my bug. This:

// zero-stuff
usidx=0;
loop(uslen,
workbuf[usidx*usfact]=srcbuf[usidx];
memset(workbuf[usidx*usfact+1],0,usfact-1);
usidx+=1;
);

... should have been this:

// zero-stuff
usidx=0;
loop(uslen,
workbuf[usidx*usfact]=srcbuf[usidx];
memset(workbuf+usidx*usfact+1,0,usfact-1);
usidx+=1;
);

The crazy thing about filters is how you can make mistakes like this and they still basically work much of the time!

I can't wait to redo all the testing I've been doing for the last two weeks with this, newly correct, implementation, to find out how well all these darn filters actually work in reality. :-)

Ouch, this was a painful one.

Thanks for all the assistance!
clepsydrae is offline   Reply With Quote
Old 12-27-2018, 08:07 AM   #12
sai'ke
Human being with feelings
 
sai'ke's Avatar
 
Join Date: Aug 2009
Location: Germoney
Posts: 195
Default

Quote:
Originally Posted by mschnell View Post
Doing more reasoning...

I am not convinced regarding the sinc (or any) filter solution.

Due to theory the samples define a unique band limited waveform and same is identical to that band limited waveform which had been sampled.
Do you mean practically or theoretically?

If frequency content is what you want to preserve, then I'm pretty certain that the full sinc is the ideal filter for rejecting the aliases.

I mean, you probably know this, but just to rehash the theory a bit. Please correct me if I'm wrong. Upsampling results in a spectrum that periodically repeats at multiples of the samplerate. The higher copies have to be rejected. Convolution (filtering) in the time domain equates to multiplication in the frequency domain and vice versa. So, one way to reject our copies is to multiply by a rectangle in the frequency domain (nulling the copies), or alternatively, convolve (filter) with a sinc (ifft of a rectangle) in the time domain. Since sincs are infinitely long, they have to be cut short. While in theory, it would be possible to just take a subset of the sinc, this typically causes a lot of problems since truncation causes distortion. Truncation in the time domain, corresponds to multiplication of the kernel by a rectangle, again leading to sharp edges. Again, sharp edge in one domain leads to ripple effects in the other. This is why the sinc is typically windowed. People typically use a Kaiser window for this, but you can also use Hamming, Blackman or some other windowing function for this.

Some people use other filters/interpolation methods, but as far as I know, this is often for performance considerations rather than theory, but please correct me if I'm wrong. If I am wrong on this, I'd be very interested in reading any references you have on filters that are more optimal and in what sense they are more optimal.

Quote:
Originally Posted by mschnell View Post
To use any kind of "normal" (electronic or software) filter, you would need to provide values for the time between two samples ("hold", linear interpolation, quadratic interpolation, spline interpolation, ...) this obviously creates a distortion to the signal to be reproduced, even if perfect band-limiting might keep the damage at bay. (See the "adding Zeros" picture in the first post here, which does not state what exactly the "interpolation filter" is supposed to be. But maybe adding zeros can be proven to cause a gain reduction as the only "distortion".)
I mean, for most of these, it's not so hard to calculate the result on the spectrum. The zero hold for instance, is a convolution with a constant block. Its frequency domain analog is again a sinc. So what you'd be doing is multiplying the frequency domain by a sinc. Not a great lowpass filter.

First order hold (linear interpolation) corresponds to upsampling then filtering with [0.5, 1, 0.5]. A triangle is nothing more than two rectangles convolved together. Convolution is multiplication in the frequency domain, so that means the triangle corresponds to a sinc squared in the frequency domain. Squaring it makes the side lobes go down faster, but it's still not a great filter.

Non-linear filtering gets ugly, because the spectral analysis toolkit only makes it possible to really study the linear properties. In these cases, I'd be quite interested in reading about what metric is even used to measure their superior properties.

Quote:
Originally Posted by mschnell View Post
I suppose there should be an algorithm that can cope with missing values. I understand that something like this is used in "oversampling D/A converters".
I am pretty certain windowed sinc or approximations thereof are what most audio resamplers employ, as S&H, linear, cubic and other interpolators tend to leave high frequency aliasing.

Cascaded comb filters are also pretty common but as far as I know, those are used mostly for performance reasons when we want to oversample by a lot (they require only addition and subtraction).

@clepsydrae
Glad you resolved the bug!

And yes, that's always the difficulty with numerical work. It's very error prone. And in the case of audio, a lot of stuff acts as a lowpass filter.

That's also why I recommended using a ground truth if you can get one.

Are you planning to release your tuner when it is finished?
__________________
[Soundcloud] | [Tracker Plugin: Thread|Github|Reapack link] | [Machine UI Plugin: Thread|Github|Reapack link] | [Filther: Github]

Last edited by sai'ke; 12-27-2018 at 08:26 AM.
sai'ke is offline   Reply With Quote
Old 12-27-2018, 09:29 PM   #13
clepsydrae
Human being with feelings
 
clepsydrae's Avatar
 
Join Date: Nov 2011
Posts: 2,236
Default

Quote:
Originally Posted by sai'ke View Post
Are you planning to release your tuner when it is finished?
Yeah... assuming it turns out to be useful. :-) I'm not yet sure if I'm prototyping a shareware thing or just a novelty.

I've noticed that if I try to upsample to 8X or 64X in a single step, the results are slightly "fuzzy" (i.e. a tiny bit of high frequency contamination), presumably because the filters are harder to generate (since the cutoff frequency is effectively 8X or 64X lower in the frequency range.)

If I upsample in stages (2X2X2X2X) it seems to go better. (This is true for my custom filters and for the upsampling library.) I still have a tiny bit of fuzz, but it's not as bad.

Does this all seem normal and correct?
clepsydrae is offline   Reply With Quote
Old 12-28-2018, 02:13 AM   #14
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Germany
Posts: 6,312
Default

I think we agree that:

We want to upsample a signal with sample rate <f> to sample rate <n*f>.

The signal is denoted by samples called s(i), i = 0… .

The the original signal is band limited to f/2 before sampling and hence the samples define a unique f/2-band-limited signal (no other f/2-band-limited signal matches all samples), which is exactly the original signal.

The (analog) <f>-restoration-hardware reproduces exactly this signal, with an imperfection that depends on the exactness of band limiting filter and the “restoration filter” provided.

An <n*f> restoration hardware needs n times as many samples to reproduce an n*f/2-band-limited signal.

The upsampled signal is denoted by samples called u(j), j = 0… , with u(n*i) hopefully close to s(i), u(j) chosen so that a perfect <n*f> restoration hardware reproduces the original <f/2>-band-limited signal (ignoring any latency introduced).

To generate the additional samples we could (e.g.):
(0) set them to zero (use no existing values, keeping u(n*i) = s(i) )
(1) do a “hold” (use one existing value, keeping u(n*i) = s(i) )
(2) do a linear interpolation (use two existing values, keeping u(n*i) = s(i) )
(3) do a quadratic interpolation (use three existing values, keeping u(n*i) = s(i) )
(4) do a cubic spline interpolation (use four existing values, keeping u(n*i) = s(i) )
(5) do a cubic spline approximation (using more than four existing values, also re-calculating the samples u(n*i) )

(“Dithering” by adding some (shaped) noise also might be discussed.)

Obviously none of this is an <n*f> representation of an f/2-band-limited signal and hence provides artifacts / distortion to the original signal in any frequency range. Hopefully the higher the id (complexity) of the algorithm, the lower the power of the distortion.

Now adding a digital filter working with the sample frequency n*f could help to reduce the distortion by reducing the power of the signal in the band between f/2 and n*f/2 (such a filter can’t do anything else). But as the signal u(i) is not a perfect <n*f> representation of a “desirable” n*f/2-band-limited signal, such a filter is prone to add artifacts below n*f/2 and even f/2, which might even have some power in a low frequency range, that the <n*f> restoration hardware reproduces as (audible) distortion.

Supposedly with a given algorithm and a given digital filter design, the amount of such distortion can be (asymptotically) calculated. (Maybe it even is zero with an infinite sinc.)


My reasoning is:

You can think the samples to be dirac pulses, each pushing the reproduction filter away from its resting value to generate a signal. With this, you don’t need a fixed sampling frequency, but samples that don’t need to “push” are obsolete and can be ignored as “missing values”. With that, the upsampling process just means do add n-1 “missing value” samples u(n*i+k) (k=1..n-1) between each defined sample u(n*i). Of course the existing (analog) <n*f>-restoration-hardware can’t cope with missing values (even though it would be possible to design such a thing in hardware) and we need to insert a digital reproduction filter, that can.

I am not fluent with digital filter design, so I just want to provide an example that will not be decently usable.

I understand that a simple first grade R/C-type (IIR-) filter is done as
v(j+1) = v(j) + a*(u(j)-v(j))
Hence a simple (Gaussian) second grade filter is done as
v(j+1) = v(j) + a*(u(j)-v(j)); w(j+1) = w(j) + a*(v(j)-w(j))

Here you see the dirac pulses u(j)-v(j) that “push” the reproduction filter.

Such a filter easily can be empowered to handle missing values:

For j = n*i →
with u(j) = s(n*i): v(j+1) = v(j) + a*(u(j)-v(j)); w(j+1) = w(j) + a*(v(j)-w(j))
else →
v(j+1) = v(j); w(j+1) = w(j) + a*(v(j)-w(j))

I gather this paradigm should easily be extensible to higher order and more steep and hence in fact usable “filters”.

Obviously, no sample of the upsampled signal w(j) will be one of the original sample values.

For me it seems like the interpolation algorithms (0).. try to create / predict such zero power dirac pulses for the missing values with increasing effort.

I suppose that such an algorithm (which I suppose is not decently called a filter) should not introduce artifacts.

What do you think ?

I suppose the CPU overhead of such a design might be rather good regarding an interpolation plus standard filter design with similar predicted distortion values.

Do you think there is a way to calculate the distortion spectrum of such a thing ?

-Michael (being very astonished if that is a new idea )
__________________
www.boa-sorte.de

Last edited by mschnell; 12-28-2018 at 03:39 AM.
mschnell is online now   Reply With Quote
Old 01-02-2019, 01:38 PM   #15
sai'ke
Human being with feelings
 
sai'ke's Avatar
 
Join Date: Aug 2009
Location: Germoney
Posts: 195
Default

Quote:
Originally Posted by clepsydrae View Post
Does this all seem normal and correct?
It's hard to tell from just a description. I mean, you can look at the spectrum of a band limited signal and observe what it does?

@mschnell
To be honest, I'm having a bit of trouble following your line of reasoning (I blame lack of sleep ).

I mean, in a sense this is already done with a FIR? The missing values are zero and therefore those taps are skipped? Usually up-sampling and filtering are done in one step. So the non-zeroes are the only ones that give a push and as the filter slides over the samples, different steps are alternatingly used?
__________________
[Soundcloud] | [Tracker Plugin: Thread|Github|Reapack link] | [Machine UI Plugin: Thread|Github|Reapack link] | [Filther: Github]
sai'ke is offline   Reply With Quote
Old 01-02-2019, 03:54 PM   #16
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Germany
Posts: 6,312
Default

Maybe this can be done with FIR. (setting the bins for frequency > f/2 to zero and then doing a (reverse) FFT). "My" suggestion is based on IIR, supposedly needing a lot less CPU power.

-Michael
__________________
www.boa-sorte.de
mschnell is online now   Reply With Quote
Old 01-03-2019, 04:35 AM   #17
sai'ke
Human being with feelings
 
sai'ke's Avatar
 
Join Date: Aug 2009
Location: Germoney
Posts: 195
Default

I think there are better options for designing FIRs for this purpose. IIRs shine at low frequencies, but also distort the phase, which may or may not be a problem depending on which filter you're using.

If you go for a non-sinc interpolator due to performance constraints, there's a nice doc with some interpolation designs and their magnitude response here:
http://yehar.com/blog/wp-content/upl...09/08/deip.pdf

Actually, all this talk of upsampling/downsampling has made me want to revisit Filther's up/downsampling again. There I used a pretty crappy RBJ lowpass as up/downsampling filter. Will be fun replacing that with a decent linear phase FIR as alternative option
__________________
[Soundcloud] | [Tracker Plugin: Thread|Github|Reapack link] | [Machine UI Plugin: Thread|Github|Reapack link] | [Filther: Github]

Last edited by sai'ke; 01-03-2019 at 08:26 AM.
sai'ke is offline   Reply With Quote
Old 01-05-2019, 01:11 AM   #18
SaulT
Human being with feelings
 
Join Date: Oct 2013
Location: Seattle, WA
Posts: 748
Default

Did someone say upsampling? lol

I've had to take a break from programming for a while, but I saw the topic and had to comment. Interpolation is a topic that sits in my brain, I think about it a lot.


TLDR; you can use any lowpass filter when oversampling, including interpolators, but some are better than others and there's usually a trade-off going on. Also, don't use the "optimized" or "parabolic" interpolators from deip.pdf, they blow out your control points.


@sai'ke

I'm going to start from what I think is the beginning, bear with me, and I apologize for the wall of text, but this is what I'm passionate about in DSP.

We have a samplerate f (let's use 48 kHz and get all "real life" with this). For some reason we want to saturate a 5 kHz sine wave. Distortion adds harmonics. Specifically, at multiples of 5 kHz, so 10 kHz, 15 kHz, 20 kHz, 25 kHz, etc. But hey, wait a sec, our samplerate only goes up to 24 kHz! The 25 kHz harmonic folds over once it hits the Nyquist and appears as aliasing - a very inharmonic 23 kHz. The next harmonic, which should have been 30 kHz, becomes 18 kHz, and so on. These bear no relation to 5 kHz and will sound like poo. To keep from aliasing, we need to raise our Nyquist, and we do that by oversampling! Oversampling 48 kHz by a factor of 2 means that our Nyquist is now at 48 kHz, and those harmonics may fade away in that higher frequency space (usually higher harmonics become progressively lower in volume).

To oversample by a factor of n we zero-stuff... we insert n zeroes between the values of the original signal. So, s(0) 0 0 0 s(1) 0 0 0 s(2) 0 0 0 etc is what zero-stuffing for 4x looks like. The problem is that this actually duplicates the signal, it creates mirror images in those higher frequencies. As a result, we have to filter them out. That is why we use a lowpass filter after we zero-stuff.

To downsample we are going to throw away a bunch of values, but when we do that any frequencies higher than our original Nyquist will then be mirrored into our signal, and that's aliasing and that's bad. So we use another lowpass filter to remove as much of that higher-than-Nyquist stuff as possible. We run that lowpass filter and we throw away the values we don't need, ie we only keep every nth value. If we've done this all correctly we are back to having a bandlimited signal and everyone is happy....

....except life isn't perfect and filters aren't perfect. It becomes a tradeoff. You can use an IIR but then you have changed phase. Maybe that's an issue, maybe not. Maybe you need linear phase. So you use a FIR, except that filter performance (how steep it can be) is based on the number of taps. Use more taps and lower the amount of possible aliasing... but it also increases your latency. In a realtime application (e.g. playing through a guitar pedal emulator) you can't have FIRs with hundreds of taps, it just isn't tenable, too much latency. So you compromise. Sometimes you can use math to your advantage, though. One option is a Nyquist filter.

A Nyquist filter (aka L-th band filter) is one that puts the cutoff at f/L. Every Lth coefficient of that FIR (counting out from the middle) is zero, and we can use that to cut down on the amount of calculations we do. If you look at Tale's code (my code is an extension of his, he inspired me) you should be able to puzzle out how it works, anytime there would be a multiplication by 0 we can just ignore it.

So, we want to oversample x2. We insert a zero between each value, then we run it through an L=2 Nyquist filter. Since our Nyquist doubled by inserting those zeroes, and we've set L=2, our operating cut-off is our original Nyquist. We save a bunch of calculations, so it's cheaper to do than a full FIR, but there are disadvantages, too. I've already gone on far too long about it, but that's my summary.

Oh, but wait, interpolation. Interpolation (and it took me a long time to understand this) is basically a lowpass filter you calculate on the fly. It will filter out those mirrored frequencies I mentioned above, but not as well as a sinc-based filter. Another way to say that is that interpolation is noisy, with better algorithms being less noisy. You can see a measurement of how good the algorithms are at filtering out data, a few pages into deip.pdf.

You can actually measure this noise (or at least see pretty pictures of it) if you have Matlab. I think I have a script somewhere, if you want I'll look for it.

When it comes to algorithms in deip.pdf, if you've read this far, I have to emphasize - the parabolic 2x and all of the "optimized" algorithms create curves that DO NOT go through your control points. You must run every point through this algorithm. It's easier to simply not use them for actual interpolation. For some reason I ended up zeroing in on the 6-point 3rd-order Hermite, and I was happy with it, but I don't remember why I went with that and not B-spline.

Okay, that's quite enough for right now. My apologies for the wall of text.


Sault
SaulT is offline   Reply With Quote
Old 01-05-2019, 06:43 AM   #19
sai'ke
Human being with feelings
 
sai'ke's Avatar
 
Join Date: Aug 2009
Location: Germoney
Posts: 195
Default

@Sault;
Hey, thanks for what you wrote, looks correct to me and summarizes it quite nicely. Much better than what I wrote. But uh, why are you addressing this to me?

I was actually arguing for windowed sincs for the up-sampling in general earlier in the thread (#3 and #12), preferably via a polyphase filterbank, not polynomials. I also referred to your thread which has those LS optimized filters . Since OP wants to find peaks, I would argue that phase linearity and predictable group delay is important, so that means IIRs are off the table.

The reason I used an IIR in Filther was because of simplicity and the fact that the analog filter models themselves are IIRs already (so no linear phase for me).

@clepsydrae; Yesterday, I did implement a FIR solution for Filther as well though (in case one wants to use Filther as a graphical waveshaper only in some dynamics setting and cares about phase linearity). I was lazy and only used a Hamming window for the sinc here. Better would be a Kaiser one. I might improve this in the future. In case you want to use it, here it is:

Preparation

Run *once* when changing oversampling factor:
Code:
    sincFilterL.updateSincFilter(slider54, 8, sinc_flt,  sinc_tmp);
    sincFilterR.updateSincFilter(slider54, 8, sinc_flt2, sinc_tmp);
Sinc_flt, sinc_flt2 and sinc_tmp should be memory addresses you don't use.

Code:
// Generate windowed sinc filter at memory location FIR
// Inputs are:
//    fir   - Memory location to store windowed sinc
//    nt    - Number of taps
//    bw    - Fractional bandwidth
//     g    - Gain
function sinc(fir, nt, bw, g)
  local(a, ys, yg, yw, i, pibw2, pifc2, pidnt2, hnt)
  global()
  (
    pibw2   = 2.0*$pi*bw;
    pidnt2  = 2.0*$pi/nt;
    hnt     = 0.5*nt;
    i       = 1;
        
    loop(nt-1,
      // Sinc width
      a  = (i - hnt) * pibw2;
        
      // Sinc
      ys = (a != 0) ? sin(a)/a : 1.0;
 
      // Window gain
      yg = g * (4.0 * bw);
        
      // Hamming window (could be replaced with Kaiser in the future)
      yw = 0.54 - 0.46 * cos(i * pidnt2);
         
      // Calc FIR coeffs
      fir[i-1] = yw * yg * ys;
      
      i += 1;
    );
  );

// Generate sinc filters for a specific upsampling ratio
//
// Upsampling leads to a sample followed by N-1 zeroes. Hence 
// to compute each subsample, we only need 1/Nth of the taps.
// This is why we set up a specific filter for each subsample.
// i.e. for N=4, you get something like f1*Zn + f5*Zn-1 + ...
//
// Inputs:
//    N_in            - oversampling factor
//    tapsPerFilter   - Taps per subfilter (should be 8 in this implementation)
//    targetmem       - Location to store the coefficients
//    tmp             - Working memory
function updateSincFilter(N_in, tapsPerFilter, targetmem, tmp)
  local(nHist, iFilt, nTaps)
  instance(h0, h1, h2, h3, h4, h5, h6, coeffs, loc, N, delta)
  global()
  (
    N       = N_in;
    nHist   = tapsPerFilter;
    loc     = 0;
    coeffs  = targetmem;
    nTaps   = N*nHist;
    
    // Memory being set is conservatively large.
    memset(coeffs,0,10000);
    memset(tmp,0,10000);
    
    sinc(tmp, nTaps, .5/N, .5*N);
    
    // Divide sinc over the different filters
    iFilt = 0; // Filter idx for which subsample this filter is
    delta = 0; // Sample idx
    loop(nTaps,
      coeffs[delta + iFilt*100] = tmp[];
      iFilt += 1;
      iFilt == N ? ( iFilt = 0; delta += 1 );
      tmp += 1;
    );
  );
Update on incoming samples
Run for every incoming sample.
Code:
      sincFilterL.advanceSinc(spl0);
      sincFilterR.advanceSinc(spl1);
Code:
// Maintain input sample history. Hard-coded for speed.
function advanceSinc(sample)
  instance(h0, h1, h2, h3, h4, h5, h6, h7, coeffs, loc, N)
  global()
  local(filt, out)
  (
    h7 = h6;
    h6 = h5;
    h5 = h4;
    h4 = h3;
    h3 = h2;
    h2 = h1;
    h1 = h0;
    h0 = sample;
    loc = 0;
  );
Now you can work with subsamples as though they were coming in normally. You should call it N times (8 in this implementation) for every advanceSinc call. You can simply make a loop with getSubSample() and just loop over your subsampled signal as though it was incoming.

Code:
sincFilterL.getSubSample();
sincFilterR.getSubSample();
Code:
function getSubSample()
  instance(out, h0, h1, h2, h3, h4, h5, h6, h7, coeffs, loc, N)
  global()
  local(filt)
  (
    filt = coeffs + loc;

    out =  filt[] * h0 + filt[1] * h1 + filt[2] * h2 + filt[3] * h3 + filt[4] * h4 + filt[5] * h5 + filt[6] * h6 + filt[7] * h7;

    loc += 100;
    out
  );
Note that if the oversampling factor is an integer number here, then filter coefficient 7 is always zero, such that filt[7] and h7 can be omitted.

The group delay (in the upsampled sequence) for this thing is (.5 * 8 * N - 1).
__________________
[Soundcloud] | [Tracker Plugin: Thread|Github|Reapack link] | [Machine UI Plugin: Thread|Github|Reapack link] | [Filther: Github]

Last edited by sai'ke; 01-05-2019 at 01:03 PM.
sai'ke is offline   Reply With Quote
Old 01-06-2019, 04:15 PM   #20
clepsydrae
Human being with feelings
 
clepsydrae's Avatar
 
Join Date: Nov 2011
Posts: 2,236
Default

Thanks guys for the generous input --

Quote:
Originally Posted by sai'ke View Post
It's hard to tell from just a description. I mean, you can look at the spectrum of a band limited signal and observe what it does?
I just ran some more tests, and I must have been confused in my filter design (shocking, I know.) This time around I can upsample 8X and there's no problem.

@SaulT -- one thing I didn't quite understand was the curve at this link -- I'm interpreting that to be 44.1k upsampled to 88.2k -- if we are upsampling, don't we need to response to be negligible above the original nyquist? Both filter curves seem to be only around -5dB down at nyquist, hence my confusion.

Quote:
Originally Posted by sai'ke View Post
Since OP wants to find peaks, I would argue that phase linearity and predictable group delay is important, so that means IIRs are off the table.
Yeah, that seemed right to me, too. I settled on a FIR filter made here. I didn't bother to optimize it ("polyphase" implementation, I take it) because my application is not real-time, but I might in the future.

Thanks, Sai'ke, for the very interesting upsampling code -- it will take me some time to digest, but if/when I get around to going through it, many thanks for the time you took!
clepsydrae is offline   Reply With Quote
Old 01-08-2019, 04:49 PM   #21
SaulT
Human being with feelings
 
Join Date: Oct 2013
Location: Seattle, WA
Posts: 748
Default

Quote:
@SaulT -- one thing I didn't quite understand was the curve at this link -- I'm interpreting that to be 44.1k upsampled to 88.2k -- if we are upsampling, don't we need to response to be negligible above the original nyquist? Both filter curves seem to be only around -5dB down at nyquist, hence my confusion.
So it's all a matter of compromise. What you gain in fewer computations with Nyquist filters is that the cutoff is always -6 dB at 1/L. The more taps you add the steeper the roll-off will be above cutoff and the less roll-off you see below it, but of course the more group delay and processing power you need.

In the link you can see that at 19 points it's much steeper than 7 points but it's still going to alias. Full attenuation is 60 dB, that makes just about everything silent. People can't hear much over 15-ish kHz. So, pick a number of taps that will keep group delay and calculations low but will still keep aliasing out of audible frequencies (e.g. here 30 kHz will still alias to 14 kHz, but it will be pulled down 60 dB and that's basically silent).

The important thing to me is that it's all a trade-off. You can always add more taps, go to a higher level of oversampling, or pick a different method entirely. For instance, you can generate a FIR with a cutoff lower than Nyquist (e.g 19 or 17 kHz instead of 22.05 kHz here) and while it wouldn't give you the benefits of fewer calculations you can get a lot more attenuation because you don't have to be so steep. And of course if you aren't constrained by linear phase then you can do even more.
SaulT is offline   Reply With Quote
Old 01-08-2019, 05:20 PM   #22
clepsydrae
Human being with feelings
 
clepsydrae's Avatar
 
Join Date: Nov 2011
Posts: 2,236
Default

Ok, thanks. I didn't realize a LPF for upsampling was afforded that much wiggle room. I figured it had to happen between like 20k and nyquist and go down 120dB or something along those lines.

I've seen reference to a half-band filter... that's the same as what you describe here as an Lth-band where L==2, or a "Nyquist" filter, right?
clepsydrae is offline   Reply With Quote
Old 01-08-2019, 11:27 PM   #23
SaulT
Human being with feelings
 
Join Date: Oct 2013
Location: Seattle, WA
Posts: 748
Default

It depends on what your design requirements are. Maybe for a "mastering quality" plugin you do want that 20k to Nyquist passband with 90+ dB attenuation above Nyquist. The more relaxed your requirements the more options you'll have and less processing you'll need.

Yes, an Lth band filter where L==2 is also known as a halfband filter. The only other common term I've heard is L==4, a "quarterband" filter.

As an example, I have a quarterband filter somewhere in here where I used the design criteria where at 44.1 kHz the passband was meant to reach full attenuation by 15 kHz, since I figured most people wouldn't hear aliasing above that frequency. At 31 points, the filter gives less than a millisecond of total latency while still giving linear phase and 60 dB (IIRC) of aliasing attenuation.
SaulT is offline   Reply With Quote
Old 01-08-2019, 11:30 PM   #24
clepsydrae
Human being with feelings
 
clepsydrae's Avatar
 
Join Date: Nov 2011
Posts: 2,236
Default

Another big thanks to everyone that pitched in to help me out, here. The waters have cleared slightly. :-)
clepsydrae is offline   Reply With Quote
Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off

Forum Jump


All times are GMT -7. The time now is 03:05 AM.


Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2019, vBulletin Solutions Inc.