Analog-derived and purely digital filters?

DSP, Plugin and Host development discussion.
Post Reply New Topic
RELATED
PRODUCTS

Post

So I've noticed that a lot of digital filter designs seem to be "adaptations" of analog filters to the digital domain. Any idea why is this? Is it perhaps easier to develop filters in analog than digital?

Are there any purely digital filter designs?

Post

There is some ambiguity in the usage of the word analog. It could mean hardware analog or simply abstract continuous-time (the former kind of implies the latter but not vice versa). Abstract continuous-time filters are must easier to understand and handle than digital (that is, discrete time) ones (YMMV) and usually have natural behavior out of the box. OTOH, proper modeling of hardware analog can be close to rocket science, depending on how many analog artifacts are there.

Post

Fluky wrote:So I've noticed that a lot of digital filter designs seem to be "adaptations" of analog filters to the digital domain. Any idea why is this? Is it perhaps easier to develop filters in analog than digital?
As Vadim (z1202) said, "analog" may or may not mean it exists in hardware, but in general it's a filter in continuous time, and that's where we exist. Nature can make a filter (hollowed log, our mouth and the rest of our vocal tract. motion in water), damping in mechanical systems, and later a capacitor and resistor. Given how long digital has existed, and the fact we don't spend any time living in it, should you expect to find filters waiting for us there? :wink: Mainly, we take filters from the analog world that do the job that we know we want to do, and find a way to design the same outcome in the digital domain.
Are there any purely digital filter designs?
Take an audio stream, delay it, and add it to the original. Digital filter design. Oh, you want a specific filter, one that has the properties you want in the continuous domain? Well, you might start with the analog description (equation) of it, and transform it to the digital domain.

But of course there are purely digital design methodologies (where the line separating the two is it depends on your definitions). FIR design methods look pretty digital to me, for instance.

I think you need to read the threads here on the board, and specific papers like Vadim's work. Then you can ask specific questions about what you don't understand. With the threads you're starting up, you seem to be wanting to learn from ground zero with questions, while bypassing the large body of work that is available that would cover all that. If you wanted to learn how to build a house, you wouldn't go on a construction forum and start asking about what kind of wood to use, and how they get the stucco to stick—you'd read the construction guides, look for threads where some of your question had been answered, and be left with specific questions that the experts could help you with. No offense meant, I'm just not sure if you're aware of all of the info available. It would also be easier to help you if you said what you're after (EQ filters? Synth?).
My audio DSP blog: earlevel.com

Post

I think you need to read the threads here on the board, and specific papers like Vadim's work. Then you can ask specific questions about what you don't understand. With the threads you're starting up, you seem to be wanting to learn from ground zero with questions, while bypassing the large body of work that is available that would cover all that. If you wanted to learn how to build a house, you wouldn't go on a construction forum and start asking about what kind of wood to use, and how they get the stucco to stick—you'd read the construction guides, look for threads where some of your question had been answered, and be left with specific questions that the experts could help you with. No offense meant, I'm just not sure if you're aware of all of the info available. It would also be easier to help you if you said what you're after (EQ filters? Synth?).
This book works well for ground zero bootstrapping: http://www.dspguide.com
Vadim's papers, especially the book is not for everyone (it wasn't for me), I'd suggest looking elsewhere.
His papers are a bit more accessible than the book because they keep the math away (to some extent). Might be a good way to continue after Steven W. Smith's book.
~stratum~

Post

BTW, your blog, www.earlevel.com is also a good ground zero bootstrapper, as you tell what exists out there and why it is important. Thanks for writing that :-)
~stratum~

Post

Thanks! (I've been tied up with a weekend DSP project the past year, feel bad about the lack of new material on my blog, but will pick it up again around the end of summer.)

Yes, agree dspguide would be a good place to start. It all depends on what he's after too—if the object is a great synth filter, maybe not so helpful, but still good general background. For synth, probably Andres Simper's papers would be a better start than Vadim's book—even if the derivation is a bit much to start with, I think it's helpful to at least see the concept of recreating the topology of an analog synth-suitable filter; and some of the discussion on this board about why such a filter is better than trying to polish a direct form biquad into a synth filter...
My audio DSP blog: earlevel.com

Post

stratum wrote:Vadim's papers, especially the book is not for everyone (it wasn't for me), I'd suggest looking elsewhere.
His papers are a bit more accessible than the book because they keep the math away (to some extent).
:lol: :lol: :lol:
It's pretty ironic, because in reality it's the other way 'round. The papers assume that you already know the math and DSP in general and also can derive the implied/omitted formulas in your head. The book simply writes all that out explicitly. In fact the book is approximately equal to the first of the papers just with extra clarifications for those who do not have the respective background. :lol: :lol: :lol:

Post

It's pretty ironic, because in reality it's the other way 'round.
There is another way around: One can understand the concepts but yet may not know the math and the presentation of math does not necessarily give an understanding of the concepts.

For example, for anyone familiar with analog electronics, Andres Simper's papers are readable up to the point he starts to use Maxima/Matlab (which one it is I dont know), after that it looks like a cumbersome set of equations. For example, looking at this one http://cytomic.com/files/dsp/SvfLinearTrapOptimised.pdf there isn't much headache causing material in the first page after the introduction, after that I cannot follow the equations, but I can see that they are transfer functions for recursive filters. Admittedly Steven W. Smith is responsible for our presence here and the "weird" way we understand things, or else we wouldn't be understanding anything at all. This is similar to the way hobbists deal with analog electronics, after a while they develop an intuitive understanding and can build a stompbox or even a complex guitar amp, and by understanding the building blocks one can even create variants of existing designs, but may fail to know what a pole or zero is. There is a difference between this level of understanding of analog filters and being able to make stability analysis of them.
~stratum~

Post

Fluky wrote:Are there any purely digital filter designs?
I've done one in Excel, years ago... I had a large series of measurements and wanted to even them out a bit so some abnormal measurements don't disrupt the graph too much. But I wanted to keep away of Excel's own smoothing algo's.

So I calculated another column: value is average of twice the current plus previous plus next measurement. That's sort of a very basic & crude LFP, I always wondered how it would sound. Now to think about it, it's a very simple IR of 3 samples long with the values 0.25, 0.5 and 0.25.
We are the KVR collective. Resistance is futile. You will be assimilated. Image
My MusicCalc is served over https!!

Post

Fluky wrote:Are there any purely digital filter designs?
Yes, of course. A FIR/correlation in the general sense can be considered impractical in the analogue domain (with exceptions here and there). But even if that's too wishi washi, take a look at advanced nonlinear filters: The median filter for example, or a filter based on percentiles. They are definitely very non-"plastic and copper". ;)
Fabien from Tokyo Dawn Records

Check out my audio processors over at the Tokyo Dawn Labs!

Post

BertKoor wrote: So I calculated another column: value is average of twice the current plus previous plus next measurement. That's sort of a very basic & crude LFP, I always wondered how it would sound. Now to think about it, it's a very simple IR of 3 samples long with the values 0.25, 0.5 and 0.25.
That's a nyquist notch filter (in my own terms). It sounds, well, kinda like bandlimiting ultrasonic frequencies if you use 96khz.

Post

Here is a funky filter from long ago. Haven't found many practical uses-- Inspired by Steven W. Smith's DSPBook observation that cascaded box filters tend toward gaussian response. Along with his observation that Gaussians can be "about as good as it gets" for smoothing broadband noise while preserving transients, though the filter is not very sharp or selective.

The only time I recall using this filter in release code was a vinyl click'n'pop remover DX Plugin. It was an optional filter which the user could enable with a checkbox, to try smoothing the "frying bacon" broadband crackle noise of badly worn records. Sometimes it "seemed to have a little benefit" for that purpose.

Have tried the filter for a few other purposes, but so far (because it is smooth and non-selective) the filter seems mostly a "solution in search of a problem". :)

The idea of cascaded boxcar filters expressed "as rudimentary as it gets"-- First stage, add the current sample and the previous sample for the first boxcar 2-sample average. Second stage, add the current two-sample average (current sample and last sample) with the previous two-sample average (last sample and the sample before that). Keep going, averaging previous averages. The FPU ASM version could do up to 6 stages, because that is how many running averages could fit in the 8-value FPU stack along with needed temp vars. It was kinda conceived as a FPU ASM filter, because it ran purt fast if you set it to work processing an entire buffer in one whack. Expressed as per-sample code in high level language, not particularly fast.

The length of 1 stage is 2 samples, the total length of 2 stages is three samples, etc. I think an 8 stage cascade would only be a length of 9 samples but might be thinking about it wrong. Dunno if that is a true gaussian after many cascaded stages. Maybe it gets skewed.

I recently made an "up to 8 stage" version in js code to test it again. This ought to work in js, but I made a couple of "last minute" changes without testing in reaper, so maybe this has bugs--

Code: Select all

//advisable to only call js_malloc within js @init section
JS_MALLOC_NO_CLEARMEM = 0;
JS_MALLOC_CLEARMEM = 1;
function js_malloc(size, clearmem)
(
  js_malloc_result = next_js_malloc;
  (clearmem != 0) ?
    memset(next_js_malloc, 0.0, size);
  next_js_malloc += size;
  js_malloc_result;
);

//call in @init section
//example: o_GS.GaussianSmooth_Init(8);
function GaussianSmooth_Init(a_NStages)
(
  this.NStages = a_NStages;
  (this.NStages < 1) ?
    this.NStages = 1
  :
    (this.NStages > 8) ?
      this.NStages = 8;
  this.Avg_Ary = js_malloc(9, JS_MALLOC_CLEARMEM); 
  //InSamp goes into Avg_Ary[0], cascaded two sample averages in Avg_Ary[1] thru Avg_Ary[8]
  this.LastAvg_Ary = js_malloc(9, JS_MALLOC_CLEARMEM); 
  //Last InSamp goes into LastAvg_ary[0], last averages in higher elements
);

//call in @sample section
//example: spl0 = o_GS.GaussianSmooth_DoSamp(spl0);
function GaussianSmooth_DoSamp(a_InSamp)
(
  _GSIdx = 0;
  this.Avg_Ary[_GSIdx] = a_InSamp;
  while (_GSIdx < this.NStages)
  (
    this.Avg_Ary[_GSIdx + 1] = (this.LastAvg_Ary[_GSIdx] + this.Avg_Ary[_GSIdx]) * 0.5;
    this.LastAvg_Ary[_GSIdx] = this.Avg_Ary[_GSIdx];
    _GSIdx += 1;
  );
  _GSResult = this.Avg_ary[this.NStages]; //js returns last assignment as function result
);
Here was some testing of the js filter--

Code: Select all

Gaussian Smooth-- Cascade Eight 2-sample averages-- Low pass characteristics

=== Detail at samprate = 44100 ===
18861 Hz = -90 dB
18022 Hz = -84 dB
17455 Hz = -78 dB
16971 Hz = -72 dB
16475 Hz = -66 dB
15931 Hz = -60 dB
15341 Hz = -54 dB
14633 Hz = -48 dB
13924 Hz = -42 dB
13074 Hz = -36 dB
12118 Hz = -30 dB
10972 Hz = -24 dB
9667 Hz = -18 dB
7977 Hz = -12 dB
5750 Hz = -6 dB
4127 Hz = -3 dB
3302 Hz = -2 dB
2316 Hz = -1 dB
1568 Hz = -0.5 dB
1386 Hz = -0.4 dB
1191 Hz = -0.3 dB
934 Hz = -0.2 dB
542 Hz = -0.1 dB

=== Generic any samprate ===
About -90 dB at Fs / 2
About -48 dB at Fs / 3
About -26.6 dB at Fs / 4
About -6 dB at Fs / 8
About -1.5 dB at Fs / 16
About -1 dB at Fs / 19
About -0.1 dB at Fs / 81

=====================

Gaussian Smooth -- Cascade Four 2-sample averages

3144 Hz = -1 dB
749 HZ = -0.1 dB

About -90 dB at Fs / 2
About -24 dB at Fs / 3
About -14 dB at Fs / 4
About -6 dB at Fs / 5.5
About -3 dB at Fs / 8
About -1 dB at Fs / 14
About -0.1 dB at Fs / 59

====================

Gaussian Smooth -- Cascade Three 2-sample averages

About -90 dB at Fs / 2
About -18 dB at Fs / 3
About -9 dB at Fs / 4
About -2.1 dB at Fs / 8
About -1 dB at Fs / 12


======================

Gaussian Smooth -- Cascade Two 2-sample averages

About -90 dB at Fs / 2
About -12 dB at Fs / 3
About -6 dB at Fs / 4
About -1.8 dB at Fs / 8
About -1 dB at Fs / 9.5

=======================

Smooth One 2-sample average

About -57 dB at Fs / 2
About -6 dB at Fs / 3
About -3 dB at Fs / 4
About -0.7 dB at Fs / 8
And the old FPU ASM delphi object--

Code: Select all

TGaussianBoxFilt = class(TObject) //Possibly useful fast filter for high-freq denoising
  private
  protected
  public
  LastVals: array [0..5] of double;
  NumStages: integer;
  Mult: double;

  constructor Create(NewNumStages: integer); //pass 1 thru 6
  procedure BeforeDestruction; override;

  procedure Reset(NewNumStages: integer; ResetOnly: boolean);
  function ProcessSample(input: double): double;
  procedure ProcessMonoBuffer(InFloatPtr: PFloatBuf; InputNumBytes: integer);
  procedure ProcessStereoBufferLeft(InFloatPtr: PFloatBuf; InputNumBytes: integer);
  procedure ProcessStereoBufferRight(InFloatPtr: PFloatBuf; InputNumBytes: integer);
end;

begin
  inherited Create;
  Reset(NewNumStages, false);
end;

procedure TGaussianBoxFilt.BeforeDestruction;
begin
  //
end;

procedure TGaussianBoxFilt.Reset(NewNumStages: integer; ResetOnly: boolean);
var
  i: integer;
begin
  for i := 0 to 5 do
    LastVals[i] := 0.0;
  if not ResetOnly then
  begin
    if NewNumStages < 1 then
      NewNumStages := 1
    else if NewNumStages > 6 then
      NewNumStages := 6;
    NumStages := NewNumStages;
    Mult := 1.0/(1 shl NumStages);
  end;
end;

function TGaussianBoxFilt.ProcessSample(input: double): double;
var
  TmpOut: double;
  i: integer;
begin
  TmpOut := input;
  for i := 0 to NumStages - 1 do
  begin
    TmpOut := input + LastVals[i];
    LastVals[i] := input;
    input := TmpOut;
  end;
  result := TmpOut * Mult;
end;

procedure TGaussianBoxFilt.ProcessMonoBuffer(InFloatPtr: PFloatBuf; InputNumBytes: integer);
var
  SelfPtr: pointer;
  PIndx: integer;
  PInTop: integer;
begin
  asm
    mov     SelfPtr,eax         //save pointer to self, eax holds self on entry
  end;
  PIndx := integer(InFloatPtr); //use integer for pointers, fewer typecasts
  PInTop := PIndx + InputNumBytes;
  asm
    push    esi
    push    edi
    push    ebx

    mov     esi,PIndx
    mov     edi,PInTop
    mov     ebx,SelfPtr
    mov     edx,[ebx].NumStages  //edx holds NumStages for compare

    fldz                         //st0 = dummy input val
    lea     eax,[ebx].LastVals
    fld     double([eax])        //st0 = LastVals[0], st1 = input
    add     eax,8
    fld     double([eax])        //st0 = LastVals[1], st1 = LastVals[0], st2 = input
    add     eax,8
    fld     double([eax])        //st0 = LastVals[2], st1 = LastVals[1], st2 = LastVals[0], st3 = input
    add     eax,8
    fld     double([eax])        //st0 = LastVals[3], st1 = LastVals[2], st2 = LastVals[1], st3 = LastVals[0], st4 = input
    add     eax,8
    fld     double([eax])        //st0 = LastVals[4], st1 = LastVals[3], st2 = LastVals[2], st3 = LastVals[1], st4 = LastVals[0], st5 = input
    add     eax,8
    fld     double([eax])        //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0], st6 = input
    jmp     @StartLoop

@ContinueLoop:
    //input := PSingle(PIndx)^;
    fld     single([esi])        //st0 = input, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = input
    fst     st(7)                //st0 = input, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = newinput
    //TmpOut := input + LastVals[0];
    fadd    st(0),st(6)          //st0 = TmpOut = input + LastVals[0], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = newinput
    //LastVals[0] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(6)                //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0] (updated), st6 = newinput = TmpOut
    cmp     edx,1                //see if NumStages = 1
    jz      @WriteOutVal

    //TmpOut := input + LastVals[1];
    fld     st(6)                //st0 = newinput = OldTmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    fadd    st(0),st(5)          //st0 = TmpOut = newinput + LastVals[1], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    //LastVals[1] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(5)                //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1] (updated), st5 = LastVals[0], st6 = input = TmpOut
    cmp     edx,2                //see if NumStages = 2
    jz      @WriteOutVal

    //TmpOut := input + LastVals[2];
    fld     st(6)                //st0 = newinput = OldTmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    fadd    st(0),st(4)          //st0 = TmpOut = newinput + LastVals[2], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    //LastVals[2] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(4)                //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3], st3 = LastVals[2] (updated), st4 = LastVals[1], st5 = LastVals[0], st6 = input = TmpOut
    cmp     edx,3                //see if NumStages = 3
    jz      @WriteOutVal

    //TmpOut := input + LastVals[3];
    fld     st(6)                //st0 = newinput = OldTmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    fadd    st(0),st(3)          //st0 = TmpOut = newinput + LastVals[3], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    //LastVals[3] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(3)                //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3] (updated), st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0], st6 = input = TmpOut
    cmp     edx,4                //see if NumStages = 3
    jz      @WriteOutVal

    //TmpOut := input + LastVals[4];
    fld     st(6)                //st0 = newinput = OldTmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    fadd    st(0),st(2)          //st0 = TmpOut = newinput + LastVals[4], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    //LastVals[4] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(2)                //st0 = LastVals[5], st1 = LastVals[4] (updated), st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0], st6 = input = TmpOut
    cmp     edx,5                //see if NumStages = 3
    jz      @WriteOutVal

    //TmpOut := input + LastVals[5];
    fld     st(6)                //st0 = newinput = OldTmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    fadd    st(0),st(1)          //st0 = TmpOut = newinput + LastVals[5], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    //LastVals[5] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(1)                //st0 = LastVals[5], st1 = LastVals[4] (updated), st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0], st6 = input = TmpOut

@WriteOutVal:
    //PSingle(PIndx)^ := TmpOut * Mult;
    fld     st(6)                //st0 = TmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = input = TmpOut
    fmul    [ebx].Mult         //st0 = TmpOut * Mult, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = input = TmpOut
    fstp    single([esi])        //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0], st6 = oldinput
    add     esi,4                //inc esi

@StartLoop:
    cmp     esi,edi              //check if done
    jb      @ContinueLoop        //if done, drop to ExitCleanup

@ExitCleanup:
    fstp    double[eax]         //st0 = LastVals[4], st1 = LastVals[3], st2 = LastVals[2], st3 = LastVals[1], st4 = LastVals[0], st5 = oldinput
    sub     eax,8
    fstp    double[eax]         //st0 = LastVals[3], st1 = LastVals[2], st2 = LastVals[1], st3 = LastVals[0], st4 = oldinput
    sub     eax,8
    fstp    double[eax]         //st0 = LastVals[2], st1 = LastVals[1], st2 = LastVals[0], st3 = oldinput
    sub     eax,8
    fstp    double[eax]         //st0 = LastVals[1], st1 = LastVals[0], st2 = oldinput
    sub     eax,8
    fstp    double[eax]         //st0 = LastVals[0], st1 = oldinput
    sub     eax,8
    fstp    double[eax]         //st0 = oldinput
    fstp    st(0)               //FPU cleaned up

    pop     ebx
    pop     edi
    pop     esi
  end;
  {
  //pas code version
  while PIndx < PInTop do
  begin
    input := PSingle(PIndx)^;
    for i := 0 to NumStages - 1 do
    begin
      TmpOut := input + LastVals[i];
      LastVals[i] := input;
      input := TmpOut;
    end;
    PSingle(PIndx)^ := TmpOut * Mult;
    PIndx := PIndx + sizeof(single);
  end;
  }
end;

procedure TGaussianBoxFilt.ProcessStereoBufferLeft(InFloatPtr: PFloatBuf; InputNumBytes: integer);
var
  SelfPtr: pointer;
  PIndx: integer;
  PInTop: integer;
begin
  asm
    mov     SelfPtr,eax         //save pointer to self, eax holds self on entry
  end;
  PIndx := integer(InFloatPtr); //use integer for pointers, fewer typecasts
  PInTop := PIndx + InputNumBytes;
  asm
    push    esi
    push    edi
    push    ebx

    mov     esi,PIndx
    mov     edi,PInTop
    mov     ebx,SelfPtr
    mov     edx,[ebx].NumStages  //edx holds NumStages for compare

    fldz                         //st0 = dummy input val
    lea     eax,[ebx].LastVals
    fld     double([eax])        //st0 = LastVals[0], st1 = input
    add     eax,8
    fld     double([eax])        //st0 = LastVals[1], st1 = LastVals[0], st2 = input
    add     eax,8
    fld     double([eax])        //st0 = LastVals[2], st1 = LastVals[1], st2 = LastVals[0], st3 = input
    add     eax,8
    fld     double([eax])        //st0 = LastVals[3], st1 = LastVals[2], st2 = LastVals[1], st3 = LastVals[0], st4 = input
    add     eax,8
    fld     double([eax])        //st0 = LastVals[4], st1 = LastVals[3], st2 = LastVals[2], st3 = LastVals[1], st4 = LastVals[0], st5 = input
    add     eax,8
    fld     double([eax])        //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0], st6 = input
    jmp     @StartLoop

@ContinueLoop:
    //input := PSingle(PIndx)^;
    fld     single([esi])        //st0 = input, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = input
    fst     st(7)                //st0 = input, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = newinput
    //TmpOut := input + LastVals[0];
    fadd    st(0),st(6)          //st0 = TmpOut = input + LastVals[0], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = newinput
    //LastVals[0] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(6)                //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0] (updated), st6 = newinput = TmpOut
    cmp     edx,1                //see if NumStages = 1
    jz      @WriteOutVal

    //TmpOut := input + LastVals[1];
    fld     st(6)                //st0 = newinput = OldTmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    fadd    st(0),st(5)          //st0 = TmpOut = newinput + LastVals[1], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    //LastVals[1] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(5)                //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1] (updated), st5 = LastVals[0], st6 = input = TmpOut
    cmp     edx,2                //see if NumStages = 2
    jz      @WriteOutVal

    //TmpOut := input + LastVals[2];
    fld     st(6)                //st0 = newinput = OldTmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    fadd    st(0),st(4)          //st0 = TmpOut = newinput + LastVals[2], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    //LastVals[2] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(4)                //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3], st3 = LastVals[2] (updated), st4 = LastVals[1], st5 = LastVals[0], st6 = input = TmpOut
    cmp     edx,3                //see if NumStages = 3
    jz      @WriteOutVal

    //TmpOut := input + LastVals[3];
    fld     st(6)                //st0 = newinput = OldTmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    fadd    st(0),st(3)          //st0 = TmpOut = newinput + LastVals[3], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    //LastVals[3] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(3)                //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3] (updated), st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0], st6 = input = TmpOut
    cmp     edx,4                //see if NumStages = 3
    jz      @WriteOutVal

    //TmpOut := input + LastVals[4];
    fld     st(6)                //st0 = newinput = OldTmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    fadd    st(0),st(2)          //st0 = TmpOut = newinput + LastVals[4], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    //LastVals[4] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(2)                //st0 = LastVals[5], st1 = LastVals[4] (updated), st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0], st6 = input = TmpOut
    cmp     edx,5                //see if NumStages = 3
    jz      @WriteOutVal

    //TmpOut := input + LastVals[5];
    fld     st(6)                //st0 = newinput = OldTmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    fadd    st(0),st(1)          //st0 = TmpOut = newinput + LastVals[5], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    //LastVals[5] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(1)                //st0 = LastVals[5], st1 = LastVals[4] (updated), st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0], st6 = input = TmpOut

@WriteOutVal:
    //PSingle(PIndx)^ := TmpOut * Mult;
    fld     st(6)                //st0 = TmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = input = TmpOut
    fmul    [ebx].Mult           //st0 = TmpOut * Mult, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = input = TmpOut
    fstp    single([esi])        //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0], st6 = oldinput
    add     esi,8                //inc esi

@StartLoop:
    cmp     esi,edi              //check if done
    jb      @ContinueLoop        //if done, drop to ExitCleanup

@ExitCleanup:
    fstp    double[eax]         //st0 = LastVals[4], st1 = LastVals[3], st2 = LastVals[2], st3 = LastVals[1], st4 = LastVals[0], st5 = oldinput
    sub     eax,8
    fstp    double[eax]         //st0 = LastVals[3], st1 = LastVals[2], st2 = LastVals[1], st3 = LastVals[0], st4 = oldinput
    sub     eax,8
    fstp    double[eax]         //st0 = LastVals[2], st1 = LastVals[1], st2 = LastVals[0], st3 = oldinput
    sub     eax,8
    fstp    double[eax]         //st0 = LastVals[1], st1 = LastVals[0], st2 = oldinput
    sub     eax,8
    fstp    double[eax]         //st0 = LastVals[0], st1 = oldinput
    sub     eax,8
    fstp    double[eax]         //st0 = oldinput
    fstp    st(0)               //FPU cleaned up

    pop     ebx
    pop     edi
    pop     esi
  end;
end;


procedure TGaussianBoxFilt.ProcessStereoBufferRight(InFloatPtr: PFloatBuf; InputNumBytes: integer);
var
  SelfPtr: pointer;
  PIndx: integer;
  PInTop: integer;
begin
  asm
    mov     SelfPtr,eax         //save pointer to self, eax holds self on entry
  end;
  PIndx := integer(InFloatPtr) + 4; //use integer for pointers, fewer typecasts
  PInTop := integer(InFloatPtr) + InputNumBytes;
  asm
    push    esi
    push    edi
    push    ebx

    mov     esi,PIndx
    mov     edi,PInTop
    mov     ebx,SelfPtr
    mov     edx,[ebx].NumStages  //edx holds NumStages for compare

    fldz                         //st0 = dummy input val
    lea     eax,[ebx].LastVals
    fld     double([eax])        //st0 = LastVals[0], st1 = input
    add     eax,8
    fld     double([eax])        //st0 = LastVals[1], st1 = LastVals[0], st2 = input
    add     eax,8
    fld     double([eax])        //st0 = LastVals[2], st1 = LastVals[1], st2 = LastVals[0], st3 = input
    add     eax,8
    fld     double([eax])        //st0 = LastVals[3], st1 = LastVals[2], st2 = LastVals[1], st3 = LastVals[0], st4 = input
    add     eax,8
    fld     double([eax])        //st0 = LastVals[4], st1 = LastVals[3], st2 = LastVals[2], st3 = LastVals[1], st4 = LastVals[0], st5 = input
    add     eax,8
    fld     double([eax])        //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0], st6 = input
    jmp     @StartLoop

@ContinueLoop:
    //input := PSingle(PIndx)^;
    fld     single([esi])        //st0 = input, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = input
    fst     st(7)                //st0 = input, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = newinput
    //TmpOut := input + LastVals[0];
    fadd    st(0),st(6)          //st0 = TmpOut = input + LastVals[0], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = newinput
    //LastVals[0] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(6)                //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0] (updated), st6 = newinput = TmpOut
    cmp     edx,1                //see if NumStages = 1
    jz      @WriteOutVal

    //TmpOut := input + LastVals[1];
    fld     st(6)                //st0 = newinput = OldTmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    fadd    st(0),st(5)          //st0 = TmpOut = newinput + LastVals[1], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    //LastVals[1] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(5)                //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1] (updated), st5 = LastVals[0], st6 = input = TmpOut
    cmp     edx,2                //see if NumStages = 2
    jz      @WriteOutVal

    //TmpOut := input + LastVals[2];
    fld     st(6)                //st0 = newinput = OldTmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    fadd    st(0),st(4)          //st0 = TmpOut = newinput + LastVals[2], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    //LastVals[2] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(4)                //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3], st3 = LastVals[2] (updated), st4 = LastVals[1], st5 = LastVals[0], st6 = input = TmpOut
    cmp     edx,3                //see if NumStages = 3
    jz      @WriteOutVal

    //TmpOut := input + LastVals[3];
    fld     st(6)                //st0 = newinput = OldTmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    fadd    st(0),st(3)          //st0 = TmpOut = newinput + LastVals[3], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    //LastVals[3] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(3)                //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3] (updated), st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0], st6 = input = TmpOut
    cmp     edx,4                //see if NumStages = 3
    jz      @WriteOutVal

    //TmpOut := input + LastVals[4];
    fld     st(6)                //st0 = newinput = OldTmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    fadd    st(0),st(2)          //st0 = TmpOut = newinput + LastVals[4], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    //LastVals[4] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(2)                //st0 = LastVals[5], st1 = LastVals[4] (updated), st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0], st6 = input = TmpOut
    cmp     edx,5                //see if NumStages = 3
    jz      @WriteOutVal

    //TmpOut := input + LastVals[5];
    fld     st(6)                //st0 = newinput = OldTmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    fadd    st(0),st(1)          //st0 = TmpOut = newinput + LastVals[5], st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = OldTmpOut
    //LastVals[5] := input;
    //input := TmpOut;
    fxch    st(7)                //st0 = newinput, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = TmpOut
    fstp    st(1)                //st0 = LastVals[5], st1 = LastVals[4] (updated), st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0], st6 = input = TmpOut

@WriteOutVal:
    //PSingle(PIndx)^ := TmpOut * Mult;
    fld     st(6)                //st0 = TmpOut, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = input = TmpOut
    fmul    [ebx].Mult         //st0 = TmpOut * Mult, st1 = LastVals[5], st2 = LastVals[4], st3 = LastVals[3], st4 = LastVals[2], st5 = LastVals[1], st6 = LastVals[0], st7 = input = TmpOut
    fstp    single([esi])        //st0 = LastVals[5], st1 = LastVals[4], st2 = LastVals[3], st3 = LastVals[2], st4 = LastVals[1], st5 = LastVals[0], st6 = oldinput
    add     esi,8                //inc esi

@StartLoop:
    cmp     esi,edi              //check if done
    jb      @ContinueLoop        //if done, drop to ExitCleanup

@ExitCleanup:
    fstp    double[eax]         //st0 = LastVals[4], st1 = LastVals[3], st2 = LastVals[2], st3 = LastVals[1], st4 = LastVals[0], st5 = oldinput
    sub     eax,8
    fstp    double[eax]         //st0 = LastVals[3], st1 = LastVals[2], st2 = LastVals[1], st3 = LastVals[0], st4 = oldinput
    sub     eax,8
    fstp    double[eax]         //st0 = LastVals[2], st1 = LastVals[1], st2 = LastVals[0], st3 = oldinput
    sub     eax,8
    fstp    double[eax]         //st0 = LastVals[1], st1 = LastVals[0], st2 = oldinput
    sub     eax,8
    fstp    double[eax]         //st0 = LastVals[0], st1 = oldinput
    sub     eax,8
    fstp    double[eax]         //st0 = oldinput
    fstp    st(0)               //FPU cleaned up

    pop     ebx
    pop     edi
    pop     esi
  end;
end;
Last edited by JCJR on Thu Jun 30, 2016 11:35 pm, edited 2 times in total.

Post

I had another app for the 'cascaded box filter' (cascaded moving average?) but never tried it.
The idea was to use it as an highpass filter (i.e. 1-lowpass(x)) for D.C. removal at high sampling rates because I wasn't sure about the behavior of even a first order IIR (being a self learner causes some unwarranted 'paranoia' sometimes ). Anyway, the idea was to convert floats to integers my multiplying them, say a large enough number to give enough amplitude resolution, then use 1-movingavg(x) as an highpass filter, and convert back to float. I had thought integers were simple enough to avoid 'strange' numeric problems that I wasn't able to understand and the moving average filter was a good fit for them.

Probably the actual problem was something else, and this is probably just another example of 'cascaded box filters' looking for a problem to solve.
~stratum~

Post Reply

Return to “DSP and Plugin Development”