|
 |
Sooooo ... I wanted to try to make an FFT filter.
Actually quite easy:
There's a couple tricks, of course. One being how the file is being
read. I read chunks 2*FFT Len... then seek back by FFT Len. This is so I
have enough extra data to do the overlap later.
The code that does the filtering isn't actually using the generated
filter (I was experimenting, and rather than mess with creating a filter
with a nice rolloff Like I set out to do, I just did a sharp cut at a
certain frequency, in this case its a high-pass filter...) The filter is
reflected from FFT_LEN to 2*FFT_LEN. This may not be necessary as I
believe the only thing it will affect is the resulting imaginary numbers
when the ifft is computed. Once the ifft is performed, the blocks are
blended using a Gaussian curve to eliminate any pops caused by abrupt
sample value changes. (This is also the reason the data is read in the
way it is, otherwise half the signal would be empty)
There is probably A LOT of improvement this could take, but its a neat
start. :D
Code follows (C#)
const int FFT_LEN = 4096;
static void Main(string[] args)
{
double[] dr = new double[FFT_LEN*2];
double[] di = new double[FFT_LEN * 2];
double[] fr = new double[FFT_LEN * 2];
FileStream fs =
System.IO.File.OpenRead(@"D:\Sample\sound.raw");
for (int i = 0; i < FFT_LEN; i++)
{
fr[i] = 1/Math.Pow(1+(i)/(double)FFT_LEN, 2);
}
int bytes = 0;
byte[] samples = new byte[FFT_LEN*2];
List<double> output = new List<double>();
int time = 0;
int ch = 0;
int bl = (int)((fs.Length/FFT_LEN) / 40) / 4;
while ((bytes = fs.Read(samples, 0, FFT_LEN*2)) > 0)
{
if(fs.Position != fs.Length)
fs.Seek(-FFT_LEN, SeekOrigin.Current);
time++;
if (time % bl == 0)
{
switch (ch)
{
case 0:
Console.Write(".");
break;
case 1:
Console.Write("\bo");
break;
case 2:
Console.Write("\bO");
break;
case 3:
Console.Write("\bX");
break;
}
ch = (ch + 1) % 4;
}
for (int i = 0; i < dr.Length; i++)
{
di[i] = 0;
if (i < bytes)
dr[i] = samples[i] - 128;
else
dr[i] = 0;
}
fft(FFT_LEN*2, dr, di, false);
for (int i = 0; i < FFT_LEN*2; i++)
{
// Convert to polar coordinates:
double m = Math.Sqrt(Math.Pow(dr[i], 2) +
Math.Pow(di[i], 2));
double p = Math.Atan2(di[i],dr[i]);
int n = i;
if (n >= FFT_LEN) n = -(FFT_LEN - (n + 1));
if (i < FFT_LEN-3000 || i > FFT_LEN+3000) m = 0;
dr[i] = m * Math.Cos(p);
di[i] = m * Math.Sin(p);
}
fft(FFT_LEN*2, dr, di, true);
if (output.Count > 0)
{
int start = output.Count - FFT_LEN;
for (int i = 0; i < FFT_LEN; i++)
{
double fn = Math.Exp(-Math.Pow(i, 2) / (2 *
Math.Pow(FFT_LEN/4, 2)));
//dr[i] += output[i + start];
dr[i] = dr[i]*(1-fn) + fn*output[i + start];
}
}
// Replace the last FFT_LEN samples with the samples in dr.
if(output.Count > 0) output.RemoveRange(output.Count -
FFT_LEN, FFT_LEN);
output.AddRange(dr);
}
fs.Close();
FileStream outF = File.OpenWrite(@"d:\sample\out.au");
BinaryWriter bw = new BinaryWriter(outF);
bw.Write(new byte[] { (byte)'.', (byte)'s', (byte)'n',
(byte)'d' });
bw.Write(BE(24));
bw.Write(BE(output.Count));
bw.Write(BE(2));
bw.Write(BE(44100));
bw.Write(BE(1));
for (int i = 0; i < output.Count; i++)
{
sbyte out_byte = (sbyte)Math.Min(Math.Max(output[i],
-127), 127);
//System.Diagnostics.Debug.Print("{0} = {1}",
output[i], out_byte);
//sbyte out_byte = (sbyte)(Math.Sin(i/100.0)*127);
bw.Write(out_byte);
}
outF.Close();
//fft(16, r, i, false);
//System.Console.ReadKey();
}
static byte[] BE(int n)
{
byte[] bytes = BitConverter.GetBytes(n);
return bytes.Reverse().ToArray();
}
static void fft(int n, double[] rx, double[] ix, bool inverse)
{
if (Math.Abs(Math.Truncate(Math.Log((double)n, 2)) -
Math.Log((double)n, 2)) > 0.00001)
{
throw new ArgumentException("n must be a power of 2");
}
int hn = n / 2;
int j = 0;
// Reverse bit sort
// Essentially, the indices bit order is reversed.
for (int i = 0; i < n - 1; i++)
{
if (j > i)
{
double tr = rx[j];
double ti = ix[j];
rx[j] = rx[i];
ix[j] = ix[i];
rx[i] = tr;
ix[i] = ti;
}
int k = hn;
while (k <= j)
{
j -= k;
k >>= 1;
}
j += k;
}
// FFT
//
int m = (int)Math.Log(n, 2);
int l1 = 1;
// Each stage in the recomposition
for (int l = 0; l < m; l++)
{
int l2 = l1;
l1 <<= 1;
double ur = 1;
double ui = 0;
double sr = (double)Math.Cos(Math.PI / l2);
double si = (double)Math.Sin(Math.PI / l2);
if (inverse) si = -si;
double tr;
double ti;
// Each Value in that state
for (j = 0; j < l2; j++)
{
for (int i = j; i < n; i += l1)
{
int i1 = i + l2;
tr = rx[i1] * ur - ix[i1] * ui;
ti = rx[i1] * ui + ix[i1] * ur;
rx[i1] = rx[i] - tr;
ix[i1] = ix[i] - ti;
rx[i] = rx[i] + tr;
ix[i] = ix[i] + ti;
}
tr = ur;
ur = tr * sr - ui * si;
ui = tr * si + ui * sr;
}
}
if (inverse)
{
for (int i = 0; i < n; i++)
{
rx[i] /= n;
ix[i] /= n;
}
}
}
--
~Mike
Post a reply to this message
|
|
|
|
Mike Raiford wrote:
> Sooooo ... I wanted to try to make an FFT filter.
You're aware that you can perform linear filtering operations without
using an FFT at all, right?
An FFT convolution is faster when the filter kernel is large, however.
> Actually quite easy:
>
> There is probably A LOT of improvement this could take, but its a neat
> start. :D
I sat here and thought about posting the 8 lines or so of Haskell it
would take to do the same task. And then I realised that Warp would
complain that it's not processing the data in-place, so it has different
performance characteristics, and then I'd have to build both your
version and my version and do detailed timing statistics, which would
require me to install the MS .NET runtime and some sort of C#
compiler... Nah. I have better things to do with my Friday afternoon. ;-)
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
|
|
On 8/27/2010 12:02 PM, Orchid XP v8 wrote:
> Mike Raiford wrote:
>> Sooooo ... I wanted to try to make an FFT filter.
>
> You're aware that you can perform linear filtering operations without
> using an FFT at all, right?
>
Well, yeah... of course. There's dozens of polynomials out there to do
it, too. But try designing an arbitrary filter with that. :D
The kicker is, all of those filter functions can be built using
capacitors, resistors and/or inductors.
> An FFT convolution is faster when the filter kernel is large, however.
Indeed. But using FFT makes filter design easier.
... And using FFT also makes doing a vocoder rather simple.
> I sat here and thought about posting the 8 lines or so of Haskell it
> would take to do the same task. And then I realised that Warp would
> complain that it's not processing the data in-place, so it has different
> performance characteristics, and then I'd have to build both your
> version and my version and do detailed timing statistics, which would
> require me to install the MS .NET runtime and some sort of C#
> compiler... Nah. I have better things to do with my Friday afternoon. ;-)
Thbbbbbbt.
--
~Mike
Post a reply to this message
|
|
null> wrote: