Implement a near brick-wall filter for non-realtime data

I'm not a DSP programmer, just an ordinary C++/C# guy. I'd like to
figure a simple way to process a .WAV file for doing very steep
low-pass and high-pass functions. Hopefully very little ripple within
passband. No realtime constraints, and the frequencies would be fixed,
so it probably amounts to implementation of a stock algorithm. There's
no actual spec for slope, as I'm just experimenting.  Just 'steep as
possible' for reasonable coding effort.

I'll read up on whatever is necessary, but it's difficult to tell
where to start. It's probably a very basic problem to most here. Any
leads?  Maybe there's a site that will compute coefficients?

On Tue, 01 Apr 2014 23:11:01 -0400, Bob wrote:

> I'm not a DSP programmer, just an ordinary C++/C# guy. I'd like to
> figure a simple way to process a .WAV file for doing very steep low-pass
> and high-pass functions. Hopefully very little ripple within passband.
> No realtime constraints, and the frequencies would be fixed,
> so it probably amounts to implementation of a stock algorithm. There's
> no actual spec for slope, as I'm just experimenting.  Just 'steep as
> possible' for reasonable coding effort.
> 
> I'll read up on whatever is necessary, but it's difficult to tell where
> to start. It's probably a very basic problem to most here. Any leads? 
> Maybe there's a site that will compute coefficients?

In these degenerate times, if the file is short enough you just read the 
whole damn thing into a computer, do an FFT, window, then do an IFFT.

If it's longer than that, then doing some sort of overlapping FFT (which 
I've never had to do, because I don't work with files that long!) would 
probably be easiest.

It'll sound terrible when you're done, of course.

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

On Tue, 01 Apr 2014 22:20:52 -0500, Tim Wescott
<tim@seemywebsite.really> wrote:

>On Tue, 01 Apr 2014 23:11:01 -0400, Bob wrote:
>
>> I'm not a DSP programmer, just an ordinary C++/C# guy. I'd like to
>> figure a simple way to process a .WAV file for doing very steep low-pass
>> and high-pass functions. Hopefully very little ripple within passband.
>> No realtime constraints, and the frequencies would be fixed,
>> so it probably amounts to implementation of a stock algorithm. There's
>> no actual spec for slope, as I'm just experimenting.  Just 'steep as
>> possible' for reasonable coding effort.
>> 
>> I'll read up on whatever is necessary, but it's difficult to tell where
>> to start. It's probably a very basic problem to most here. Any leads? 
>> Maybe there's a site that will compute coefficients?
>
>In these degenerate times, if the file is short enough you just read the 
>whole damn thing into a computer, do an FFT, window, then do an IFFT.
>
>If it's longer than that, then doing some sort of overlapping FFT (which 
>I've never had to do, because I don't work with files that long!) would 
>probably be easiest.
>
>It'll sound terrible when you're done, of course.

I was liking your reply up until that last bit, Tim.

Even at 24 bit / 192k sample rate, a 3 minute file should take 'only'
about 207MB in RAM, so that approach seems do-able. But I'm not clear
on how that would work phase-wise (otherwise FFT would seem to be a
great data compression method for music).  I do like the idea of
infinite attenuation, with no phase shift, if that is practical. (But
if so, why are filters still used for audio apps?)

So...why would it sound terrible, aside from losing part of the
frequency spectrum?  Or are you referring to a side-effect of the
round trip through  FFT/IFFT.

On 2.4.14 08:44, Bob wrote:
> On Tue, 01 Apr 2014 22:20:52 -0500, Tim Wescott
> <tim@seemywebsite.really> wrote:
>
>> On Tue, 01 Apr 2014 23:11:01 -0400, Bob wrote:
>>
>>> I'm not a DSP programmer, just an ordinary C++/C# guy. I'd like to
>>> figure a simple way to process a .WAV file for doing very steep low-pass
>>> and high-pass functions. Hopefully very little ripple within passband.
>>> No realtime constraints, and the frequencies would be fixed,
>>> so it probably amounts to implementation of a stock algorithm. There's
>>> no actual spec for slope, as I'm just experimenting.  Just 'steep as
>>> possible' for reasonable coding effort.
>>>
>>> I'll read up on whatever is necessary, but it's difficult to tell where
>>> to start. It's probably a very basic problem to most here. Any leads?
>>> Maybe there's a site that will compute coefficients?
>>
>> In these degenerate times, if the file is short enough you just read the
>> whole damn thing into a computer, do an FFT, window, then do an IFFT.
>>
>> If it's longer than that, then doing some sort of overlapping FFT (which
>> I've never had to do, because I don't work with files that long!) would
>> probably be easiest.
>>
>> It'll sound terrible when you're done, of course.
>
> I was liking your reply up until that last bit, Tim.
>
> Even at 24 bit / 192k sample rate, a 3 minute file should take 'only'
> about 207MB in RAM, so that approach seems do-able. But I'm not clear
> on how that would work phase-wise (otherwise FFT would seem to be a
> great data compression method for music).  I do like the idea of
> infinite attenuation, with no phase shift, if that is practical. (But
> if so, why are filters still used for audio apps?)
>
> So...why would it sound terrible, aside from losing part of the
> frequency spectrum?  Or are you referring to a side-effect of the
> round trip through  FFT/IFFT.


You'll start hearing the superposition of the missing harmonics
as ringing. To demonstrate, start with the well-known Fourier
series of a perfect square wave and look at the sums of the
lower harmonics only. The resultants are quite un-square.

Another demonstration (though not quite the same as a filter)
would be to put the speaker at the end of a length of water
pipe and listening at the other end.

< nag mode >
A different story is the 24/192 audio, see:
<http://people.xiph.org/~xiphmont/demo/neil-young.html#toc_1ch>.
</ nag mode off >

For preserving the time relationships in the filtered signal,
consider a constant-delay (symmetric, linear phase) FIR filter.
Please note that a constant-delay filter does not change the
phase of different frequencies by the same amount.

-- 

Tauno Voipio

Bob wrote:
> I'm not a DSP programmer, just an ordinary C++/C# guy. I'd like to
> figure a simple way to process a .WAV file for doing very steep
> low-pass and high-pass functions. Hopefully very little ripple within
> passband. No realtime constraints, and the frequencies would be fixed,
> so it probably amounts to implementation of a stock algorithm. There's
> no actual spec for slope, as I'm just experimenting.  Just 'steep as
> possible' for reasonable coding effort.
>
> I'll read up on whatever is necessary, but it's difficult to tell
> where to start. It's probably a very basic problem to most here. Any
> leads?  Maybe there's a site that will compute coefficients?
>

http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html

--
Les Cargill

Bob wrote:
> On Tue, 01 Apr 2014 22:20:52 -0500, Tim Wescott
> <tim@seemywebsite.really> wrote:
>
>> On Tue, 01 Apr 2014 23:11:01 -0400, Bob wrote:
>>
>>> I'm not a DSP programmer, just an ordinary C++/C# guy. I'd like to
>>> figure a simple way to process a .WAV file for doing very steep low-pass
>>> and high-pass functions. Hopefully very little ripple within passband.
>>> No realtime constraints, and the frequencies would be fixed,
>>> so it probably amounts to implementation of a stock algorithm. There's
>>> no actual spec for slope, as I'm just experimenting.  Just 'steep as
>>> possible' for reasonable coding effort.
>>>
>>> I'll read up on whatever is necessary, but it's difficult to tell where
>>> to start. It's probably a very basic problem to most here. Any leads?
>>> Maybe there's a site that will compute coefficients?
>>
>> In these degenerate times, if the file is short enough you just read the
>> whole damn thing into a computer, do an FFT, window, then do an IFFT.
>>
>> If it's longer than that, then doing some sort of overlapping FFT (which
>> I've never had to do, because I don't work with files that long!) would
>> probably be easiest.
>>
>> It'll sound terrible when you're done, of course.
>
> I was liking your reply up until that last bit, Tim.
>
> Even at 24 bit / 192k sample rate, a 3 minute file should take 'only'
> about 207MB in RAM, so that approach seems do-able. But I'm not clear
> on how that would work phase-wise (otherwise FFT would seem to be a
> great data compression method for music).  I do like the idea of
> infinite attenuation, with no phase shift, if that is practical. (But
> if so, why are filters still used for audio apps?)
>

There are a lot of different ways to make filters for digital audio.

> So...why would it sound terrible,

Because you're "bucket stomping" the FFT data. The discontinuity in the
FFT version of the data will be expressed as artifacts  when
you reverse FFT it.



> aside from losing part of the
> frequency spectrum?  Or are you referring to a side-effect of the
> round trip through  FFT/IFFT.
>

Newp. Roundtripping should not be bad with FFTW using 64 bit float
samples.

I will stop short of saying "it's lossless" but it's not far off
from lossless if it is lossy. But perhaps I misremember.

-- 
Les Cargill

On Tue, 01 Apr 2014 23:11:01 -0400, Bob <bob@nope.com>
wrote:

>I'm not a DSP programmer, just an ordinary C++/C# guy. I'd like to
>figure a simple way to process a .WAV file for doing very steep
>low-pass and high-pass functions. Hopefully very little ripple within
>passband. No realtime constraints, and the frequencies would be fixed,
>so it probably amounts to implementation of a stock algorithm. There's
>no actual spec for slope, as I'm just experimenting.  Just 'steep as
>possible' for reasonable coding effort.
>
>I'll read up on whatever is necessary, but it's difficult to tell
>where to start. It's probably a very basic problem to most here. Any
>leads?  Maybe there's a site that will compute coefficients?


 I'd vote for a simple FIR filter.  If you're not familar
with one of these, it is just a "tapped delay line"
consisting of a long chain of 1-sample delays, with the
output of each stage multiplied by a particular factor for
that stage, and all summed together to get a single output
sample.  Then the next input sample goes in, pushes
everything up one position, and repeat.

"Simple" here may be misleading, because it can have
hundreds of taps, but the design is simple. The tap weights
are essentially the impulse response of the filter you want,
which in your case is trivial because it's a rectangular
function... you don't need to do an (I)FFT to get it.  

That would make a true "brick wall" filter, but that's not
what you want because it would exhibit "Gibbs phenomenon"
.. apparent "overshoot" and "ringing".  (Those are
descriptive terms from resonant-type filters that are
applied here because the response sorta looks the same, if
you don't notice that the FIR filter seems to also
anticipate the next transition.)    

The cause of Gibbs phenomenon is using a finite set of taps,
but using more doesn't help (perhaps surprisingly).  But the
cure is pretty simple:  Instead of using the set of computed
tap weights as-is, you apply a window function to taper the
ends gently to zero... exactly the same as when windowing an
FFT, and you can use the same functions.   There are lots of
options, but I'd recommend Hann.

I use this method to generate band-limited noise from
broadband noise (random number generator).  I use up to 512
taps, which are really only 256 unique weights mirrored
about the center.  You can use arrays to index the weights
and data points into a multiplier that feeds an accumulator.
The delay line can be done with a circular buffer, so you
don't really need to move any data along, just keep
overwriting the oldest point with the newest, and advance
the pointer.

The advantage of this over an FFT is that it can handle any
length of data easily.   A single FFT for an entire WAV file
isn't generally practical, and even though you can in
principle use multiple smaller FFTs with
windowed-overlap-and-add, I think the FIR is simpler.

Best regards,





Bob Masta
 
              DAQARTA  v7.50
   Data AcQuisition And Real-Time Analysis
              www.daqarta.com
Scope, Spectrum, Spectrogram, Sound Level Meter
 Frequency Counter, Pitch Track, Pitch-to-MIDI 
   FREE Signal Generator, DaqMusiq generator    
          Science with your sound card!

On Wed, 02 Apr 2014 01:44:41 -0500, Les Cargill wrote:

> Bob wrote:
>> On Tue, 01 Apr 2014 22:20:52 -0500, Tim Wescott
>> <tim@seemywebsite.really> wrote:
>>
>>> On Tue, 01 Apr 2014 23:11:01 -0400, Bob wrote:
>>>
>>>> I'm not a DSP programmer, just an ordinary C++/C# guy. I'd like to
>>>> figure a simple way to process a .WAV file for doing very steep
>>>> low-pass and high-pass functions. Hopefully very little ripple within
>>>> passband. No realtime constraints, and the frequencies would be
>>>> fixed,
>>>> so it probably amounts to implementation of a stock algorithm.
>>>> There's no actual spec for slope, as I'm just experimenting.  Just
>>>> 'steep as possible' for reasonable coding effort.
>>>>
>>>> I'll read up on whatever is necessary, but it's difficult to tell
>>>> where to start. It's probably a very basic problem to most here. Any
>>>> leads? Maybe there's a site that will compute coefficients?
>>>
>>> In these degenerate times, if the file is short enough you just read
>>> the whole damn thing into a computer, do an FFT, window, then do an
>>> IFFT.
>>>
>>> If it's longer than that, then doing some sort of overlapping FFT
>>> (which I've never had to do, because I don't work with files that
>>> long!) would probably be easiest.
>>>
>>> It'll sound terrible when you're done, of course.
>>
>> I was liking your reply up until that last bit, Tim.
>>
>> Even at 24 bit / 192k sample rate, a 3 minute file should take 'only'
>> about 207MB in RAM, so that approach seems do-able. But I'm not clear
>> on how that would work phase-wise (otherwise FFT would seem to be a
>> great data compression method for music).  I do like the idea of
>> infinite attenuation, with no phase shift, if that is practical. (But
>> if so, why are filters still used for audio apps?)
>>
>>
> There are a lot of different ways to make filters for digital audio.
> 
>> So...why would it sound terrible,
> 
> Because you're "bucket stomping" the FFT data. The discontinuity in the
> FFT version of the data will be expressed as artifacts  when you reverse
> FFT it.
> 
> 
> 
>> aside from losing part of the frequency spectrum?  Or are you referring
>> to a side-effect of the round trip through  FFT/IFFT.
>>
>>
> Newp. Roundtripping should not be bad with FFTW using 64 bit float
> samples.
> 
> I will stop short of saying "it's lossless" but it's not far off from
> lossless if it is lossy. But perhaps I misremember.

Given infinite precision the FFT is lossless.  Since the IFFT is 
basically the FFT with one or two multiplication factors changed, it's 
lossless, too.

-- 
Tim Wescott
Control system and signal processing consulting
www.wescottdesign.com

On Wed, 02 Apr 2014 01:44:10 -0400, Bob wrote:

> On Tue, 01 Apr 2014 22:20:52 -0500, Tim Wescott
> <tim@seemywebsite.really> wrote:
> 
>>On Tue, 01 Apr 2014 23:11:01 -0400, Bob wrote:
>>
>>> I'm not a DSP programmer, just an ordinary C++/C# guy. I'd like to
>>> figure a simple way to process a .WAV file for doing very steep
>>> low-pass and high-pass functions. Hopefully very little ripple within
>>> passband. No realtime constraints, and the frequencies would be fixed,
>>> so it probably amounts to implementation of a stock algorithm. There's
>>> no actual spec for slope, as I'm just experimenting.  Just 'steep as
>>> possible' for reasonable coding effort.
>>> 
>>> I'll read up on whatever is necessary, but it's difficult to tell
>>> where to start. It's probably a very basic problem to most here. Any
>>> leads? Maybe there's a site that will compute coefficients?
>>
>>In these degenerate times, if the file is short enough you just read the
>>whole damn thing into a computer, do an FFT, window, then do an IFFT.
>>
>>If it's longer than that, then doing some sort of overlapping FFT (which
>>I've never had to do, because I don't work with files that long!) would
>>probably be easiest.
>>
>>It'll sound terrible when you're done, of course.
> 
> I was liking your reply

I aim to please

> up until that last bit, Tim.

but I'm not God.

> 
> Even at 24 bit / 192k sample rate, a 3 minute file should take 'only'
> about 207MB in RAM, so that approach seems do-able. But I'm not clear on
> how that would work phase-wise (otherwise FFT would seem to be a great
> data compression method for music).  I do like the idea of infinite
> attenuation, with no phase shift, if that is practical. (But if so, why
> are filters still used for audio apps?)
> 
> So...why would it sound terrible, aside from losing part of the
> frequency spectrum?  Or are you referring to a side-effect of the round
> trip through  FFT/IFFT.

First, you need to try it and listen to it.  If you just "stomp 
buckets" (thanks Les) then it'll cause ringing.  You can reduce, but not 
eliminate, the ringing with a shaped window (a flat top with smooth curvy 
bits down to the floor works well).

You may like the result, in which case you're home free.

You'll be implementing the "symmetric FIR" filter that Tauno mentioned -- 
basically, the FFT/IFFT pair will closely emulate any filtering that you 
can dream up.

As to why filters are still used in audio apps -- I'm not an audio guy, 
but I do know that sometimes those filters are implemented by doing FFT/
IFFT pairs in blocks.  There's a whole ocean of knowledge about how you 
do this in a continuous fashion so the joints between blocks isn't 
audible, and I've never even stuck my toe in there.  So I can't tell you 
exactly how to do it.

-- 
Tim Wescott
Control system and signal processing consulting
www.wescottdesign.com

On Wed, 02 Apr 2014 13:37:41 GMT, N0Spam@daqarta.com (Bob
Masta) wrote:

>"Simple" here may be misleading, because it can have
>hundreds of taps, but the design is simple. The tap weights
>are essentially the impulse response of the filter you want,
>which in your case is trivial because it's a rectangular
>function... you don't need to do an (I)FFT to get it.  

Here's some working code, in my Daqarta macro language, but
hopefully readable.  This is macro code I wrote to emulate
what I was doing in assembler, to insure I got the identical
results for this simple rectangular filter before
experimenting with complex shapes.  So 'L.0.BandRise' is the
rising corner frequency (ie. the low-pass cutoff) in Hz for
the Left stream 0 noise band generator.  The '65536 * 2' is
a scale factor to match the integer-based real-time scheme.
Buf1[] is a buffer that holds all the tap weights. 

The first FIR is the style you see most commonly, with and
odd number of taps (an even number, half above and half
below zero, which is the odd tap). Here I'm using 255-0-255,
for 511 total.  Typically, FIR taps are numbered with 0 at
the center, so here I've just shifted the index values up so
0 is at index 256.

R=L.0.BandRise
F=L.0.BandFall
W=(F-R) * 2 / SmplRate
Buf1[256]=W * 65536 * 2             ;Center tap weight
R=R * 2 * pi / SmplRate             ;Step, radians/sample
F=F * 2 * pi / SmplRate
A=R                                ;Accums start on 1st step
B=F
UN=1                            ;255-0-255 taps (511 total)
WHILE.UN=<256
    X=(sin(B) - sin(A)) / (UN * pi) * 65536 * 2
    Buf1[256+UN]=X                  ;Upper taps
    Buf1[256-UN]=X                  ;Mirrored lower taps
    A=A+R                           ;Next accum steps
    B=B+F
    UN=UN+1                         ;Next tap
WEND.

Here's a variant with an even number of taps (512), still
centered about 0, which is in fact missing:

R=L.0.BandRise * 2 * pi / SmplRate     ;Step,radians/sample
F=L.0.BandFall * 2 * pi / SmplRate
A=R/2                             ;Accums start at 0.5 step
B=F/2
UN=0                                    ;Tap counter
WHILE.UN=<256                ;Taps per side (512 total)
    X=(sin(B) - sin(A)) / ((UN + 0.5) * pi) * 65536 * 2
    Buf1[256+UN]=X                      ;Upper taps
    Buf1[255-UN]=X                      ;Mirrored lower taps
    A=A+R                               ;Next accum steps
    B=B+F
    UN=UN+1                             ;Next tap
WEND.

Note that this doesn't show the window function which is
applied to these later.  The Hann window which I recommended
is just a raised cosine function, such that it is unity at
the 0 tap and tapers to zero at the end taps. 

Note that all this is done ahead of time and stored, so the
actual filtering just pulls the values from the array and
multiplies by the passing data points.

Hope this helps!



Bob Masta
 
              DAQARTA  v7.50
   Data AcQuisition And Real-Time Analysis
              www.daqarta.com
Scope, Spectrum, Spectrogram, Sound Level Meter
 Frequency Counter, Pitch Track, Pitch-to-MIDI 
   FREE Signal Generator, DaqMusiq generator    
          Science with your sound card!