Analog-derived and purely digital filters?
-
- KVRian
- Topic Starter
- 1097 posts since 28 May, 2010 from Finland
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?
Are there any purely digital filter designs?
-
- KVRAF
- 1607 posts since 12 Apr, 2002
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.
-
- KVRian
- 653 posts since 4 Apr, 2010
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? 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.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?
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.Are there any purely digital filter designs?
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
-
- KVRAF
- 2256 posts since 29 May, 2012
This book works well for ground zero bootstrapping: http://www.dspguide.comI 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?).
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~
-
- KVRAF
- 2256 posts since 29 May, 2012
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~
-
- KVRian
- 653 posts since 4 Apr, 2010
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...
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
-
- KVRAF
- 1607 posts since 12 Apr, 2002
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).
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.
-
- KVRAF
- 2256 posts since 29 May, 2012
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.It's pretty ironic, because in reality it's the other way 'round.
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~
- KVRAF
- 15279 posts since 8 Mar, 2005 from Utrecht, Holland
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.Fluky wrote:Are there any purely digital filter designs?
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.
My MusicCalc is served over https!!
My MusicCalc is served over https!!
- KVRian
- 1169 posts since 24 Feb, 2012
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".Fluky wrote:Are there any purely digital filter designs?
Fabien from Tokyo Dawn Records
Check out my audio processors over at the Tokyo Dawn Labs!
Check out my audio processors over at the Tokyo Dawn Labs!
-
- KVRAF
- 7402 posts since 17 Feb, 2005
That's a nyquist notch filter (in my own terms). It sounds, well, kinda like bandlimiting ultrasonic frequencies if you use 96khz.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.
-
- KVRAF
- 3080 posts since 17 Apr, 2005 from S.E. TN
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--
Here was some testing of the js filter--
And the old FPU ASM delphi object--
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
);
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
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.
-
- KVRAF
- 2256 posts since 29 May, 2012
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.
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~