|
 |
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
|
|
