Easiest way to generate minimum phase upsampling FIR filters?

[Kraku]

I want to generate a minimum phase 4x upsampling filter.
What's the easiest way to do that?

If it was just a regular linear phase FIR, that would be a quick and simple task. But I have zero experience creating minimum phase filters.

One idea I came up with was to simply send an impulse through a bunch of IIR low pass filters and record the output into a FIR and then use that.

Are there better options, as we're basically talking about a very specific use case here, so there could be ways to generate such FIR kernels easily from scratch without having to go deep into all the "how to convert linear phase filter into a minimum phase filter" topic.

[mystran]

One idea I came up with was to simply send an impulse through a bunch of IIR low pass filters and record the output into a FIR and then use that.

Works, but it's fairly inefficient.

Are there better options, as we're basically talking about a very specific use case here, so there could be ways to generate such FIR kernels easily from scratch without having to go deep into all the "how to convert linear phase filter into a minimum phase filter" topic.

The easiest approach to converting a linear-phase FIR into minimum-phase is take FFT (padding it to a few times the original length is a good idea), take logarithm of magnitude (use some tiny bias to keep it from going to -inf), take IFFT of the result (log cepstrum), zero negative time indexes, multiply positive indexes by two, take FFT, take complex exponent and finally take IFFT back to time domain.

It's not perfect, but with large enough FFT size it gives perfectly usable results.

[diroxe7660]

I can't speak for SciPy's minimum phase converter specifically, but I've had good results with its FIR designers - maybe worth a look.

https://docs.scipy.org/doc/scipy/refere ... phase.html

[mystran]

I can't speak for SciPy's minimum phase converter specifically, but I've had good results with its FIR designers - maybe worth a look.

https://docs.scipy.org/doc/scipy/refere ... phase.html

As far as I can guess from the documentation, those are probably variations of the same basic log-cepstrum approach I outlined above.

A few variations are possible, I'm a bit hazy on details, because I've mostly just memorized the one I know works "well enough" but if I recall correctly, when we start with a linear-phase filter, we will have pairs of zeroes on both sides of the unit circle and zeroing out one side of the log-cepstrum gets rid of the zeroes outside (or inside, depending on what you zero, but it doesn't really matter 'cos you can just time-reverse the final kernel).

If we just drop half the zeroes though, the final minimum-phase response will have only half the attenuation of the original linear-phase kernel. You can design for this, but I don't particularly like this approach, because twice the attenuation in the original kernel means taking logarithms of smaller numbers (ie. the logarithms will be further away from zero) which is numerically worse for FFT precision purposes.

In theory, the solution is to reflect one side of the cepstrum to the other side (and therefore reflect the zeroes across the unit circle) which preserves the response. I'm not sure if I ever got truly good results this way though (might be me screwing up, might be numerical issues) and if we're starting with log-magnitude cepstrum taking the magnitude essentially forces the original FIR to be linear-phase no matter what and we're guaranteed symmetry in the cepstrum, so just zeroing one half and multiplying the other half by 2 (leave the index zero alone) does the job with less opportunities to get something wrong.

In practice, expect the final filter to be slightly different (eg. stop-band notches tend to move a bit), but if the conversion is done correctly (and with enough zero-padding) then the differences should be fairly negligible for practical purposes.

[martinvicanek]

Why not implement the minimum phase filter as a recursive filter? IIR filters can be polyphase and multirate for best efficiency. I have obtained good results following this resource:
R.A. Valenzuela and A.G. Constantinides,
Digital signal processing schemes forefficient interpolation and decimation
https://www.academia.edu/36770195/Digit ... decimation
They also provide a very effcient recipe for calculating the filter coefficients.

[mystran]

Why not implement the minimum phase filter as a recursive filter? IIR filters can be polyphase and multirate for best efficiency. I have obtained good results following this resource:
R.A. Valenzuela and A.G. Constantinides,
Digital signal processing schemes forefficient interpolation and decimation
https://www.academia.edu/36770195/Digit ... decimation
They also provide a very effcient recipe for calculating the filter coefficients.

Half-band (2x) IIRs are indeed quite efficient and used to be very popular, though loop-carried dependencies can make them less attractive these days at least on desktop platforms (ie. FIRs are easy to SIMD and the only loop dependency is the accumulators and every sample is independent) and the efficiency falls at higher oversampling rates compared to a polyphase FIR.

There are situations where I might choose these, there are situations where I might even brute-force an IIR, but ... a whole lot of time a polyphase FIR is straight-forward and efficient.

[FilterEverything]

Minimum phase systems have their magnitude and phase responses linked together. Essentially for every magnitude response, there is one and only one correct (minimum) phase response.
Here are the equations describing the relationship between them:
https://en.wikipedia.org/wiki/Minimum_p ... e_response

So it is possible to calculate the particular phase response that should go with the particular magnitude response, as long as you can program the calculations from the wiki page.

As I understand it, all you have to do is take the natural logarithm of the magnitude response first, apply the Hilbert transform to that, and then multiply the result by -1.

Looking at the wiki page of the Hilbert transform, it seems like it would be a good idea to take the double sided FFT of the magnitude response (so an FFT of an FFT?), calculate the phases from that, shift those phases by +90degrees for the negative frequencies and -90degrees for the positive frequencies and IFFT the result back.

Of course I have never done any of this, I don't actually know what the hell a Hilbert transform is, all I know is just some LTI system theory.

[Kraku]

I was hoping there might be some simple formula like Sinc(x), but a minimum phase version.

[joestash]

AUTO-ADMIN: Non-MP3, WAV, OGG, SoundCloud, YouTube, Vimeo, Twitter and Facebook links in this post have been protected automatically. Once the member reaches 5 posts the links will function as normal.

I'm a bit late to the party, but I've done quite a few minimum-phase resampling implementations.

If you want a polyphase FIR, the Hilbert algorithm described here works well: https://users.ece.utexas.edu/~bevans/pa ... nPhase.pdf (https://users.ece.utexas.edu/~bevans/papers/2000/minPhase/minPhase.pdf)

For example, here is a minimum-phase 64-tap per phase upsample-by-4, with +-0.01dB ripple and -108dB stopband

64tap_upby4.png

And here are the coefficients:

Code: Select all (#)

	 2.20740021e-04f,  1.31632867e-03f,  4.84611101e-03f,  1.37169155e-02f, 
	 3.24296783e-02f,  6.67732732e-02f,  1.22664755e-01f,  2.04048815e-01f, 
	 3.10170218e-01f,  4.33010696e-01f,  5.56035640e-01f,  6.55393008e-01f, 
	 7.04188436e-01f,  6.79453182e-01f,  5.70194401e-01f,  3.83938703e-01f, 
	 1.48960871e-01f, -8.97438592e-02f, -2.80796705e-01f, -3.80754833e-01f, 
	-3.68095174e-01f, -2.51854169e-01f, -7.10758348e-02f,  1.15901741e-01f, 
	 2.49586826e-01f,  2.89026942e-01f,  2.25972708e-01f,  8.79006642e-02f, 
	-7.18631643e-02f, -1.93777439e-01f, -2.34317998e-01f, -1.82257209e-01f, 
	-6.24538830e-02f,  7.45731751e-02f,  1.73269627e-01f,  1.95554283e-01f, 
	 1.35953602e-01f,  2.30549519e-02f, -9.33130460e-02f, -1.63680152e-01f, 
	-1.60016550e-01f, -8.72014606e-02f,  2.00186850e-02f,  1.13084977e-01f, 
	 1.51437066e-01f,  1.20298451e-01f,  3.65321687e-02f, -5.98815625e-02f, 
	-1.24640160e-01f, -1.29435348e-01f, -7.44496742e-02f,  1.25772170e-02f, 
	 9.00165171e-02f,  1.22135043e-01f,  9.56215629e-02f,  2.51642730e-02f, 
	-5.43140157e-02f, -1.04794421e-01f, -1.03321302e-01f, -5.24988556e-02f, 
	 2.15224242e-02f,  8.23112249e-02f,  1.01016926e-01f,  7.00858641e-02f, 
	 6.28379523e-03f, -5.81846677e-02f, -9.18387772e-02f, -7.92868946e-02f, 
	-2.82853601e-02f,  3.47672712e-02f,  7.83888528e-02f,  8.17028412e-02f, 
	 4.44437993e-02f, -1.35323975e-02f, -6.27165245e-02f, -7.89209994e-02f, 
	-5.51768340e-02f, -4.69228730e-03f,  4.63677010e-02f,  7.23878699e-02f, 
	 6.11436033e-02f,  1.95378301e-02f, -3.04601319e-02f, -6.33547077e-02f, 
	-6.31065727e-02f, -3.09510084e-02f,  1.57620452e-02f,  5.28636139e-02f, 
	 6.18453607e-02f,  3.90880191e-02f, -2.76416304e-03f, -4.17550320e-02f, 
	-5.81049869e-02f, -4.42368837e-02f, -8.25840727e-03f,  3.06855614e-02f, 
	 5.25665709e-02f,  4.67615681e-02f,  1.71947917e-02f, -2.01498563e-02f, 
	-4.58324422e-02f, -4.70621913e-02f, -2.40565387e-02f,  1.05032311e-02f, 
	 3.84203804e-02f,  4.55470076e-02f,  2.89439769e-02f, -1.98327882e-03f, 
	-3.07635329e-02f, -4.26129611e-02f, -3.20193206e-02f, -5.27026538e-03f, 
	 2.32138000e-02f,  3.86324768e-02f,  3.34852390e-02f,  1.11976625e-02f, 
	-1.60472702e-02f, -3.39448046e-02f, -3.35678390e-02f, -1.58033061e-02f, 
	 9.47082528e-03f,  2.88507239e-02f,  3.25031908e-02f,  1.91417280e-02f, 
	-3.62934743e-03f, -2.36097458e-02f, -3.05267446e-02f, -2.13053148e-02f, 
	-1.38682084e-03f,  1.84392069e-02f,  2.78651148e-02f,  2.24135717e-02f, 
	 5.53439200e-03f, -1.35148066e-02f, -2.47298415e-02f, -2.26037953e-02f, 
	-8.80978356e-03f,  8.97226228e-03f,  2.13128213e-02f,  2.20230355e-02f, 
	 1.12420105e-02f, -4.90982165e-03f, -1.77831715e-02f, -2.08212550e-02f, 
	-1.28859150e-02f,  1.39143336e-03f,  1.42853279e-02f,  1.91456076e-02f, 
	 1.38158015e-02f,  1.54960380e-03f, -1.09382097e-02f, -1.71357724e-02f, 
	-1.41195566e-02f, -3.90666091e-03f,  7.83531088e-03f,  1.49202822e-02f, 
	 1.38932945e-02f,  5.69561076e-03f, -5.04556831e-03f, -1.26137878e-02f, 
	-1.32365958e-02f, -6.95059654e-03f,  2.61488626e-03f,  1.03151954e-02f, 
	 1.22483657e-02f,  7.71985576e-03f, -5.68180489e-04f, -8.10660852e-03f, 
	-1.10233462e-02f, -8.06170559e-03f, -1.08817746e-03f,  6.05299506e-03f, 
	 9.64928996e-03f,  8.04078160e-03f,  2.36364832e-03f, -4.20249824e-03f, 
	-8.20479846e-03f, -7.72461504e-03f, -3.28065002e-03f,  2.58729383e-03f, 
	 6.75780049e-03f,  7.18060722e-03f,  3.87164501e-03f, -1.22489923e-03f, 
	-5.36464528e-03f, -6.47345424e-03f, -4.17628543e-03f,  1.19833189e-04f, 
	 4.06975339e-03f,  5.66304736e-03f,  4.23870132e-03f,  7.34479608e-04f, 
	-2.90577890e-03f, -4.80287070e-03f, -4.10502791e-03f, -1.35366237e-03f, 
	 1.89418260e-03f,  3.93886651e-03f,  3.82119590e-03f,  1.76026477e-03f, 
	-1.04620552e-03f, -3.10881713e-03f, -3.43112401e-03f, -1.98176323e-03f, 
	 3.63994375e-04f,  2.34201314e-03f,  2.97509514e-03f,  2.04844031e-03f, 
	 1.57781090e-04f, -1.65969447e-03f, -2.48901113e-03f, -1.99210518e-03f, 
	-5.31062718e-04f,  1.07472453e-03f,  2.00264105e-03f,  1.84327627e-03f, 
	 7.71236820e-04f, -5.94468027e-04f, -1.54203247e-03f, -1.63356208e-03f, 
	-9.00381838e-04f,  2.15682569e-04f,  1.12178815e-03f,  1.38495793e-03f, 
	 9.33466142e-04f,  5.83647468e-05f, -7.66195462e-04f, -1.13639468e-03f, 
	-9.12424812e-04f, -2.62143676e-04f,  4.51798731e-04f,  8.69671726e-04f, 
	 8.11125065e-04f,  3.45957681e-04f, -2.64382501e-04f, -7.12554953e-04f, 
	-7.98583866e-04f, -5.16174212e-04f, -3.82295560e-05f,  3.80421362e-04f, 
	 5.28217098e-04f,  3.28928193e-04f, -1.37217235e-04f, -6.83628151e-04f, 
	-1.10670757e-03f, -1.26985215e-03f, -1.14699773e-03f, -8.14282286e-04f, 
	-4.02911929e-04f, -4.03404731e-05f,  1.93795005e-04f,  2.82699615e-04f, 
	 2.59711361e-04f,  1.80486363e-04f,  9.56467766e-05f,  3.65237480e-05f,

But a polyphase IIR implementation will likely outperform this. If you are processing stereo, you can use SIMD to compute both channels and both polyphase paths in parallel. To upsample by 4, a 2-stage implementation is best, where the first stage has a sharp cutoff, and the second stage takes advantage of the oversampling with a lower order and mild cutoff.

Here is an example of a first-stage polyphase IIR upsample-by-2 with a very sharp cutoff and -120dB of stopband:

Code: Select all (#)

class AudioInterpolate2_SSE {
    __m128 _state[7] {};

public:
    int process(const float* input, float* output, int numFrames) {

        // stopband = -120dB @ 24.9kHz (Fs=48kHz)
        const __m128 coef0 = _mm_setr_ps(0.028421895718584f, 0.028421895718584f, 0.107532521502474f, 0.107532521502474f);   // [ C0 C0 C6 C6 ]
        const __m128 coef1 = _mm_setr_ps(0.221551815160481f, 0.221551815160481f, 0.351274010835715f, 0.351274010835715f);   // [ C1 C1 C7 C7 ]
        const __m128 coef2 = _mm_setr_ps(0.480006834844170f, 0.480006834844170f, 0.596711429539968f, 0.596711429539968f);   // [ C2 C2 C8 C8 ]
        const __m128 coef3 = _mm_setr_ps(0.696224373700835f, 0.696224373700835f, 0.777842521785289f, 0.777842521785289f);   // [ C3 C3 C9 C9 ]
        const __m128 coef4 = _mm_setr_ps(0.843567928641561f, 0.843567928641561f, 0.896722332469210f, 0.896722332469210f);   // [ C4 C4 C10 C10 ]
        const __m128 coef5 = _mm_setr_ps(0.941118482600491f, 0.941118482600491f, 0.980718999134477f, 0.980718999134477f);   // [ C5 C5 C11 C11 ]

        __m128 u0 = _state[0];
        __m128 u1 = _state[1];
        __m128 u2 = _state[2];
        __m128 u3 = _state[3];
        __m128 u4 = _state[4];
        __m128 u5 = _state[5];
        __m128 u6 = _state[6];

        for (int i = 0; i < numFrames; i++) {

            // path [ 0, 0, 1, 1 ]
            __m128 u = _mm_castpd_ps(_mm_load_sd((const double *)&input[NCH * i]));
            u = _mm_movelh_ps(u, u);    // [ l0 r0 l0 r0 ]

            u0 = _mm_fmadd_ps(coef0, _mm_sub_ps(u, u1), u0);
            u1 = _mm_fmadd_ps(coef1, _mm_sub_ps(u0, u2), u1);
            u2 = _mm_fmadd_ps(coef2, _mm_sub_ps(u1, u3), u2);
            u3 = _mm_fmadd_ps(coef3, _mm_sub_ps(u2, u4), u3);
            u4 = _mm_fmadd_ps(coef4, _mm_sub_ps(u3, u5), u4);
            u5 = _mm_fmadd_ps(coef5, _mm_sub_ps(u4, u6), u5);

            u6 = u5;
            u5 = u4;
            u4 = u3;
            u3 = u2;
            u2 = u1;
            u1 = u0;
            u0 = u;

            _mm_storeu_ps(&output[NCH * 2*i], u6);
        }

        _state[0] = u0;
        _state[1] = u1;
        _state[2] = u2;
        _state[3] = u3;
        _state[4] = u4;
        _state[5] = u5;
        _state[6] = u6;

        return 2 * numFrames;
    }
};