# Wave filters

15 replies to this topic

### #1 chaz13

chaz13

GMC Member

• GMC Member
• 3827 posts
• Version:Unknown

Posted 25 April 2012 - 09:29 PM

Hey.

I'm working on a sort of acoustic modelling program (and have a synth/sampler package I'm working on which this applies to as well). For starters, here's some basic info as to what I'm doing.

I'll be creating/changing sounds stored in WAV files. These files are basically a header, with a bunch of data that follows. This data is what creates the wave shape; each byte stores a snapshot of the height of the wave at a certain time, so the long list of data basically just stores the shape of the sound wave, which is directly describes the motion of the speaker cone.

Now, the frequency of a sound is defined by the number of waves a second, Hz. Any general sound will have a wide range of frequencies simultaneously, some harmonics (the same musical note in octaves, that is 12 semitones higher) and other noise.

This is where my problem comes in. I want to go about programming a filter, which is basically an algorithm that removes/accentuates certain frequencies in the waveform- For example, cutting all frequencies above or below a certain point. ALL sound manipulation programs have this feature, even ones as simple as audacity.

My issue is that I can't understand how you go from a number of pieces of data that explain the height of a wave, i.e the amplitude in time, to removing certain frequencies of noise from the waveform. My first thought is to simply subtract sine waves (the fundamental wave shape that has NO other noise or harmonics, only contains the sound at the exact frequency of the wave) above/below the cutoff point, for each possibly frequency.

E.g, if I wanted a High Pass filter at 30Hz (i.e, any sound below 30Hz is removed), I'd subtract sine waves with a frequency of 29Hz to 1Hz. however, I find it hard to believe this is a valid method as it'd be slow and there'd be a large amount of noise left over from intermittent frequencies (20.5Hz, or whatever).

The main problem is I'm having an extremely hard time finding any documentation on how to go about programming such a filter, which I find odd as it's a widely used function (not only in sound applications, in dealing with any kind of wave forms).

Hopefully someone here will have an idea on how to go about doing this, or possibly know of some decent documentation. in fact, failing that, hopefully someone would have the slightest idea of what I've just said.

EDIT: After reading around a little, it seems I need to transfer the data from the time to frequency domain. (So amplitude over frequency, rather than time). Now I just need to work out how to do that- I think it might require some kind of snapshot like system but I'm not entirely sure how.

Edited by chaz13, 25 April 2012 - 09:38 PM.

• 0

### #2 IceMetalPunk

IceMetalPunk

InfiniteIMPerfection

• Retired Staff
• 9260 posts
• Version:Unknown

Posted 25 April 2012 - 09:44 PM

You need to use a Fourier transform to do the conversion. And that's about all I know on that subject . Sadly, I never did learn about them in any more detail than what they're used for (which is odd, since as a computer science major, at least 3 courses I've already taken were supposed to cover them and didn't). Hopefully, since you have motivation for learning about them, the Wikipedia article will steer you in the right direction.

This may also be more applicable to your situation: Fourier series. It breaks the overall wave down into its individual sine and cosine waves, with the arguments of the trig functions of course indicating the frequencies. Once you do that, you can simply subtract the terms with frequencies outside your desired range and recompose the remaining terms.

-IMP

Edited by IceMetalPunk, 25 April 2012 - 09:47 PM.
More specific information found

• 0

### #3 loverock125

loverock125

GMC Member

• GMC Member
• 1606 posts
• Version:GM8

Posted 25 April 2012 - 09:50 PM

I can't think of any other way other than subtracting sine waves at the specific frequency.

I'm just curious though, how are you going to find all the frequencies the sound uses, in Game Maker? By comparing the amplitude at each point?

Edit: I'm not sure if you have access to this site, but if you have, then it may be proven useful:
http://staffweb.cms.gre.ac.uk/~sp02/introductiontosignals/lecture.html

Later it was shown that any non periodic waveform could also be described as the sum of sine waves if all frequencies (not just harmonically related ones were included.)

(which is odd, since as a computer science major, at least 3 courses I've already taken were supposed to cover them and didn't)

Don't only first year students of computer science study a course related to waves? I thought this module had to do more with hardware and programming than signals.

Edited by loverock125, 25 April 2012 - 09:56 PM.

• 0

### #4 IceMetalPunk

IceMetalPunk

InfiniteIMPerfection

• Retired Staff
• 9260 posts
• Version:Unknown

Posted 25 April 2012 - 10:18 PM

I'm just curious though, how are you going to find all the frequencies the sound uses, in Game Maker? By comparing the amplitude at each point?

A .wav file contains information about the sample rate, which is samples/sec. Samples_read/sample_rate = time, so you can count the number of peaks (local maxima of the samples) per unit time and voila, there's your frequency. Of course, decomposing a full sound into its individual sinusoidal waves first is more complicated, but that's what Fourier series are for, I'm told.

(which is odd, since as a computer science major, at least 3 courses I've already taken were supposed to cover them and didn't)

Don't only first year students of computer science study a course related to waves? I thought this module had to do more with hardware and programming than signals.

I was supposed to learn it in math class freshman year as a prerequisite for other classes. I never did. Then I was supposed to learn it in one of my comp sci classes, and never did. Now I'm taking an intro to algorithms class, and although we don't actually need to know Fourier series for the course, the course description lists them as "may be useful to know for this class". Haven't needed them yet, though I'm sure they can be very useful .

-IMP
• 0

### #5 chaz13

chaz13

GMC Member

• GMC Member
• 3827 posts
• Version:Unknown

Posted 25 April 2012 - 10:28 PM

Thanks for the help. It's awesome you know what I'm talking about, hopefully you didn't find my crude explanations patronizing I'm starting an acoustical engineering degree in October so I'm doing a little fiddling of my own to get it kick started.

I managed to find the Fourier transform page just before you linked me to it, but your brief explanation has already given me more of an understanding than I've gleaned from surfing the web. It seems to be all old engineering/physics pages from the 90's, electronics/hardware implementations, or Wikipedia which makes the simplest ideas far too overcomplicated (or perhaps I'm just a little dense, who knows? )

I'll try to find more information on Fourier transformations and how to apply them, then give it a shot I suppose. If anyone has any more ideas/documentation I'd be very grateful for the help.
• 0

### #6 IceMetalPunk

IceMetalPunk

InfiniteIMPerfection

• Retired Staff
• 9260 posts
• Version:Unknown

Posted 25 April 2012 - 10:39 PM

Well, there's always the Simple Wiki version . It has a source link that (I think) isn't in the original Wikipedia footnotes, so that might help.

-IMP
• 0

### #7 chaz13

chaz13

GMC Member

• GMC Member
• 3827 posts
• Version:Unknown

Posted 25 April 2012 - 11:11 PM

I don't suppose there's a normal-wiki-without-all-the-subject-specific-vocabulary?

From what I understand so far, Fourier transformations take snapshots of the waveform, and on each snapshot uses a kind of expansion of a waveform called a Fourier series (based on euler's identity) which fairly accurately represents the original waveform in terms of sine and cosine functions (where the coefficient of each trig function is related to the frequency). And this can then be used to represent the wave in terms of frequency against amplitude. This can then be reversed using the same function?

Still not sure how to apply the transformation but I'm getting there.

Edited by chaz13, 25 April 2012 - 11:20 PM.

• 0

### #8 Desert Dog

Desert Dog

GMC Member

• Global Moderators
• 6409 posts
• Version:Unknown

Posted 25 April 2012 - 11:53 PM

Ha, this is way out of my league. I could barely read a .wav file.

A bit of googling (aimed at stackoverflow) lead me down to this, though:
http://www.dspguide.com/

Which looks pretty extreme, but you can read the chapters online, and it does have a section on Fourier transformation (Chapter 10 )

http://www.dspdimens...g-using-the-ft/

I'm rusty as(I NEED that patronizing explanation! ), but I believe the Pitch is the same as the frequency (the x part of the sine wave) correct?

Other than this, I just want to say good luck, I'll watch the topic with interest. I doubt I'll be able to assist much.
• 0

### #9 IceMetalPunk

IceMetalPunk

InfiniteIMPerfection

• Retired Staff
• 9260 posts
• Version:Unknown

Posted 26 April 2012 - 12:05 AM

Ha, this is way out of my league. I could barely read a .wav file. :tongue:

A .wav file is ridiculously simple to read. Once you get past the headers, it's literally just "the byte value is how far the speaker moves forward now".

I'm rusty as(I NEED that patronizing explanation! :tongue: ), but I believe the Pitch is the same as the frequency (the x part of the sine wave) correct?

Pitch and frequency are basically the same, yeah, but frequency is not the horizontal component of the wave . It's the number of times the wave peaks in 1 second. So if it peaks 5 thousand times in one second, it's a 5kHz wave. Common sampling rates for WAV are upwards of 4 thousand samples per second, so it's easy to record frequencies in the entire range of human hearing (200Hz to 20kHz).

-IMP
• 0

### #10 Desert Dog

Desert Dog

GMC Member

• Global Moderators
• 6409 posts
• Version:Unknown

Posted 26 April 2012 - 01:34 AM

A .wav file is ridiculously simple to read. Once you get past the headers, it's literally just "the byte value is how far the speaker moves forward now".

I know, I know. I was referring to an earlier topic of chaz's, about a year ago, where we messed about, and eventually sorted out how to write (and read) .wav files.
• 0

### #11 xot

xot

media multimixer

• Global Moderators
• 4650 posts
• Version:GM:Studio

Posted 26 April 2012 - 02:21 PM

Nice DSP link, DD. That should definitely come in handy for my own audio projects.
• 0

### #12 chaz13

chaz13

GMC Member

• GMC Member
• 3827 posts
• Version:Unknown

Posted 26 April 2012 - 10:22 PM

That is a fantastic link. Just bought a couple of awesome engineering textbooks too which apparently have fairly large sections on DSP which should be handy. I'm still working on fully understanding the process of using the Fourier transformation but I'm slowly getting somewhere.
• 0

### #13 Yourself

Yourself

The Ultimate Pronoun

• Retired Staff
• 7341 posts
• Version:Unknown

Posted 28 April 2012 - 06:11 AM

I haven't looked at the book link posted above, but I have dealt with digital signal processing recently and I can share how I understand it's approached (and thus spare you the task of trying to piece it together through Wikipedia and bits of Matlab documentation like I did). One possibility is to actually take the Fourier transform as others have suggested. This is only viable if you have all of the data available from the start and the Fourier transform allows you to apply a filter essentially on all the data at once. If you're streaming data, I don't know that that can be made to work. This is where digital filters come in. And lots of math.

We usually operate in the Laplace domain for continuous filters with something called the transfer function. In the most basic sense, our filter has an input, x(t), that's a function of time, and an output, y(t) that's a function of time. If we take the Laplace transform of both of these, we get functions X(s) and Y(s) in the frequency domain (where s is a complex number sometimes referred to as the complex frequency). I don't want to spend too much time on what the Laplace transform is or what it actually does but it is similar to the Fourier transform (in fact, the Fourier transform is the Laplace transform evaluated for values of s along the imaginary axis, so you can think of the Laplace transform as being an expanded Fourier transform). At this point we introduce something called the transfer function, H(s), which we define as H(s) = Y(s) / X(s).

Most filters are described in terms of their transfer function. Roughly speaking, the transfer function basically tells you what the output of the filter would be if you put in a signal with a particular frequency. But don't worry about that too much, this part is really just a means to an end. The first thing we're going to start with is the transfer function of a low-pass filter. It turns out that if we have that, we can create any kind of filter we want through some simple transformations. For the time being, I'm going to assume that the cut-off frequency of our low-pass filter is 1 rad/s. Later on a warping factor will show up that will allow us to adjust this.

But, first off, the transformations from a low-pass transfer function, Hlow(s), to other transfer functions:

High pass:
$H_{high}(s) = H_{low}\left( \frac{1}{s} \right)$

Band pass:
$H_{bpass}(s) = H_{low}\left( Q ( s + \frac{1}{s} ) \right)$
Where the Q-factor is computed from the fractional bandwidth of the filter. This also assumes that the center frequency of the filter (the geometric mean of the two cutoff frequencies) is 1 (more on this later).

Band stop:
$H_{bstop}(s) = H_{low}\left( \frac{1}{Q ( s + \frac{1}{s} )} \right)$

So, what does a low-pass filter transfer function actually look like? Well, the simplest one that I know of is a Butterworth low-pass filter (this is actually a first order Butterworth low-pass filter, higher order filters are possible):

$H_{low}(s) = \frac{1}{1 + s}$

So, we can get the equivalent Butterworth high pass filter:

$H_{high}(s) = \frac{1}{1 + \frac{1}{s}} = \frac{s}{1 + s}$

I'll leave it as an exercise to do the other ones.

Okay, so now we have a transfer function, how do we take a stream of sampled numbers and filter them using this? Well, the first step is to transform this transfer function from the continuous time domain, to the discrete time domain (after all, we're working with discrete samples, the waveform is only recorded at regular intervals). For that we need something called the bilinear transform. Given some analog filter (continuous time) transfer function, Ha(s), we want to find its digital (discrete time) transfer function, Hd(z). Now, the Wikipedia article uses s = 2 / T ( z - 1 ) / ( z + 1 ) as its transform. I'm going to do something slightly different (recall when I mentioned that I was going to assume the cutoff frequency of the low-pass filter was always 1, well here's the part where we take care of that). Instead of 2 / T, I'm going to introduce a warping factor, that we'll call c for now.

As the article discusses, this transformation, since it's only an approximation, introduces a frequency warping, which causes the cut-off frequency to move. Replacing that 2 / T with c and following the algebra (easy enough since that 2 / T factor always stays in front), we arrive at the following relation:

$\omega_a = c \tan\left( \omega \frac{T}{2} \right)$

Where T is our sampling interval (the amount of time that passes between discrete samples in our input). Now, this relation tells us that, for some frequency, ω, the frequency that actually corresponds to in the analog filter is ωa. Our frequencies get warped! Now, what we'd like is for some specific frequency to map to another specific frequency. For example, since our analog low-pass filter has a cutoff frequency is 1, we'd like the desired cutoff frequency to map to 1:

$1 = c \tan\left( \omega_c \frac{T}{2} \right)$
$c = \frac{1}{\tan\left( \omega_c \frac{T}{2} \right)}$

So our high and low-pass filters have a cutoff frequency (the frequency above or below which we want to attenuate the response). But what about our band-pass and band-stop filters? Those really have two cutoff frequencies, a low and a high frequency, that defines their bands. Well, this is a good time to talk about that Q factor and the center frequency. First, the center frequency of a band-pass or band-stop filter is defined as the geometric mean of its high and low cutoff frequencies. That is:

$\omega_0 = \sqrt{\omega_1 \omega_2}$

Now, you might be wondering why it's not simply the average of the two. Well, technically it is the average of the two frequencies...in log space. If you take the logarithms of the lower frequency and the upper frequency, then the logarithm of the center frequency will be their average. Now let's get to that Q factor. The Q factor is defined as:

$Q = \frac{\omega_0}{\omega_2 - \omega_1}$

One short aside, the warping factor (that c term) for band-pass and band-stop filters will be the center frequency (as I mentioned briefly before, the transformations I used above assume that the center frequency is 1, so we want our actual center frequency to map to 1).

Okay, we still haven't full moved into the digital domain, so let's continue with our high-pass filter example. We left off at:

$H_a(s) = \frac{s}{1 + s}$

Let's apply our bilinear transform:

$H_d(z) = H_a\left(c \frac{z - 1}{z + 1}\right) = \frac{c \frac{z - 1}{z + 1}}{1 + c \frac{z - 1}{z + 1}}$

Doing some algebra:

$H_d(z) = \frac{ c z - c }{ z + 1 + c z - c } = \frac{ c z - c }{ (1 + c) z + (1 - c) } = \frac{ c - c z^{-1} }{(1 + c) + (1 - c)z^{-1}}$

Now, we have a digital transfer function. Like the continuous transfer function, it's defined as the ratio of the z-transforms (the discrete time equivalent of the Laplace transform) of the output and the input:

$H_d(z) = Y(z) / X(z) = \frac{ c - c z^{-1} }{(1 + c) + (1 - c)z^{-1}}$

We're much closer now, a little more algebra to go:

$\left[( 1 + c ) + ( 1 - c )z^{-1}\right] Y(z) = \left[ c - c z^{-1} \right] X(z)$
$(1 + c) Y(z) + (1 - c) z^{-1} Y(z) = c X(z) - c z^{-1} X(z)$

Now, we're going to take the inverse z-transform. The only tricks we need to know are that the z-transform is linear, so constant multiples (like c or 1 + c) will just pop out and that the inverse of z-n X(z) is x[k - n] where k is the current time sample (meaning that this basically refers to previous time samples). So, here we go:

$(1 + c) y[k] + (1 - c) y[k - 1] = c x[k] - c x[k - 1]$
$y[k] = \frac{1}{1 + c}\left( c x[k] - c x[k - 1] - (1 - c) y[k - 1] \right)$

So what does this tell us? This tells us that the value of the output at the current sample, k, is a combination of the current and previous inputs, x[k] and x[k - 1], and the previous output, y[k - 1]. So, to implement our high-pass filter in code, it might look like this:

var den, c;
c = 1 / tan( pi * cutoffFrequency_Hz / samplingFrequency_Hz );

output = ( c * ( input - lastInput ) - ( 1 - c ) * lastOutput ) / ( 1 + c );

lastInput = input;
lastOutput = output;


Note: I've changed the expression for c to use frequencies rather than angular frequency and interval, it results in a slightly simpler expression.

So now we can just step through the samples one at a time computing each output. One important note is that the cutoff frequency should be greater than 0 and less than half the sampling frequency. You should never try to filter frequencies above half the sampling frequency (essentially because of aliasing, discrete data simply can't capture frequencies that high). I'll leave it as an exercise for others to go through the derivations of the other filter types.

As a bit of interesting trivia, let's analyze our filter to see if it makes sense. Let's let the cutoff frequency tend towards 0. Since this is a high-pass filter, that means we're essentially telling it to not attenuate any frequencies, so our output should match our input.

Well, if cutoff frequency goes to 0, we can see that c tends towards infinity. If we take that limit in our expression for the output, we get this:

output ~= input - lastInput + lastOutput

Basically, it's taking the difference between the current and previous inputs (how much they've changed) and adding it to the previous output. This will make the output follow the input, exactly as we expect.

If, instead, we go the other direction and let c tend towards 0 (the case where we let the cutoff frequency approach half the sampling frequency), we find the following:

output ~= lastOutput

meaning that the filter doesn't ever change!
• 3

### #14 chaz13

chaz13

GMC Member

• GMC Member
• 3827 posts
• Version:Unknown

Posted 28 April 2012 - 12:41 PM

That's absolutely fantastic, thank you. I'll take a look at that in more detail and have a look at implementing it in one of my current projects to see how well I get it working!

EDIT:
It works like a dream. I'm so grateful, thanks Yourself. Now to implement all the other filter types... then combine them to make a parametric EQ!

Is there a general way of boosting certain frequencies, or would you apply a filter like the high pass above, then add the result to the original wave?

One- I'm just looking at a glance, but where you're using the frequency transformation, why does this:

$= H_a \left(j \frac{2}{T} \cdot \frac{ \left(e^{j \omega T/2} - e^{-j \omega T/2}\right) /(2j)}{\left(e^{j \omega T/2} + e^{-j \omega T/2 }\right) / 2}\right)$

Not cancel down to a hyperbolic function?

MODERATOR EDIT: Image replaced with TeX so that it is readable when using dark forum skin.

Edited by xot, 29 April 2012 - 02:59 AM.

• 0

### #15 Tepi

Tepi

GMC Member

• Global Moderators
• 4201 posts
• Version:GM8.1

Posted 28 April 2012 - 06:44 PM

why does this:

$= H_a \left(j \frac{2}{T} \cdot \frac{ \left(e^{j \omega T/2} - e^{-j \omega T/2}\right) /(2j)}{\left(e^{j \omega T/2} + e^{-j \omega T/2 }\right) / 2}\right)$

Not cancel down to a hyperbolic function?

It actually does. But it cancels down to tanh whose argument contains the imaginary unit, which is usually not practical.
• 1

### #16 Yourself

Yourself

The Ultimate Pronoun

• Retired Staff
• 7341 posts
• Version:Unknown

Posted 29 April 2012 - 02:40 AM

why does this:

$= H_a \left(j \frac{2}{T} \cdot \frac{ \left(e^{j \omega T/2} - e^{-j \omega T/2}\right) /(2j)}{\left(e^{j \omega T/2} + e^{-j \omega T/2 }\right) / 2}\right)$

Not cancel down to a hyperbolic function?

It actually does. But it cancels down to tanh whose argument contains the imaginary unit, which is usually not practical.

And tanh(i x) = i tan(x), which ends up cancelling nicely with an extra i that shows up. Here are some other relevant identities:

$\sinh(\imath x) = \imath \sin(x)$
$\cosh(\imath x) = \cos(x)$
• 1

#### 0 user(s) are reading this topic

0 members, 0 guests, 0 anonymous users