Simple DFT Lowpass

Question

I'm having some trouble making a simple lowpass filter with a DFT. In the end, I hope to be able to pitch-shift audio in real time, but as it stands I can't even get this right. I have no training in this area, I only know that FFTs change waves to frequencies and iFFTs do that back, and a couple of other things I've read. To bo honest I'm surprised it works as well as it does so far. Anyway here's the code:

        byte[] samples = new byte[20000000];
        int spos = 0;

samples is filled here with 8Bit Unsigned PCM. spos <- number of samples

        int samplesize = 128;
        int sampleCount = spos / samplesize;
        frequencies = new System.Numerics.Complex[sampleCount][];
        for (int i = 0; i < sampleCount; i++)
        {
            Console.WriteLine("Sample " + i + " / " + sampleCount);
            frequencies[i] = new System.Numerics.Complex[samplesize];
            for (int j = 0; j < samplesize; j++)
            {
                frequencies[i][j] = (float)(samples[i * samplesize + j] - 128) / 128.0f;
            }
            dft.Radix2Forward(frequencies[i], MathNet.Numerics.IntegralTransforms.FourierOptions.Default);
        }

        int shiftUp = 1000; //1khz
        int fade = 2; //8 sample fade.
        int kick = frequencies[0].Length * shiftUp / rate;

So now I've calculated a bunch of DFTs for 128 sample portions of the input. kick is (I hope) the number of samples in the DFT that span 1000Hz. I.E since frequencies.Length / 2 contains frequency amplitude data up to rate/2 Hz, then frequencies[0].Length / 2 * shiftUp / (rate / 2) = frequencies[0].Length * shiftUp / rate should give me the right value

        for (int i = 0; i < sampleCount; i++)
        {

This is the part I have trouble with. Without it, the output sounds great! This skips both index 0 and index 64. Both of these have a complex component of 0, and I recall reading somewhere that the value at index 0 was important...

            for (int j = 0; j < frequencies[i].Length; j++)
            {
                if (j == 0 || j == 64)
                    continue;
                if (j < 64)
                {
                    if (!(j < kick + 1))
                    {
                        frequencies[i][j] = 0;
                    }
                }
                else
                {
                    if (!(j - 64 > 63 - kick))
                    {
                        frequencies[i][j] = 0;
                    }
                }
            }

Finally it undoes the transform

            dft.Radix2Inverse(frequencies[i], MathNet.Numerics.IntegralTransforms.FourierOptions.Default);

...tosses it back in the samples array

            for (int j=0; j<samplesize; j++)
                samples[i * samplesize + j] = (byte)(frequencies[i][j].Real * 128.0f + 128.0f);
        }

...chucks it into a file

        BinaryWriter bw = new BinaryWriter(File.OpenWrite("sound"));
        for (int i = 0; i < spos; i++)
        {
            bw.Write(samples[i]);
        }
        bw.Close();

...then I import it into Audacity to murder my ears with artifacts.

The spectral display shows that the code works, to an extent

However there's these annoying highpitched crackling sounds that occur throughout the entire song. I've heard something about the Gibbs phenomenon and a window function, but I don't really know how to apply that here. The fade variable is my best attempt at a window function: everything past the 1000hz mark fades to 0 in 2 samples.

Any ideas?

Thanks!

Answer 1

So it turns out that I was right (yay): Every 1024 samples I was getting a click sound which made the audio sound awful. To fix this, I faded inbetween many short overlapping chunks of filtered audio. It's not fast, but it works, and I'm pretty sure this is what they mean by "windowing"

public class OggDFT
{
    int sample_length;
    byte[] samples;
    DragonOgg.MediaPlayer.OggFile f;
    int rate = 0;
    System.Numerics.Complex[][] frequencies;
    DiscreteFourierTransform dft = new DiscreteFourierTransform();
    int samplespacing = 128;
    int samplesize = 1024;
    int sampleCount;

    public void ExampleLowpass()
    {
        int shiftUp = 1000; //1khz
        int fade = 2; //8 sample fade.
        int halfsize = samplesize / 2;
        int kick = frequencies[0].Length * shiftUp / rate;
        for (int i = 0; i < sampleCount; i++)
        {
            for (int j = 0; j < frequencies[i].Length; j++)
            {
                if (j == 0 || j == halfsize)
                    continue;
                if (j < halfsize)
                {
                    if (!(j < kick + 1))
                    {
                        frequencies[i][j] = 0;
                    }
                }
                else
                {
                    if (!(j - halfsize > halfsize - 1 - kick))
                    {
                        frequencies[i][j] = 0;
                    }
                }
            }
            dft.BluesteinInverse(frequencies[i], MathNet.Numerics.IntegralTransforms.FourierOptions.Default);
        }
    }

    public OggDFT(DragonOgg.MediaPlayer.OggFile f)
    {
        Complex[] c = new Complex[10];
        for (int i = 0; i < 10; i++)
            c[i] = i;
        ShiftComplex(-2, c, 5, 10);


        this.f = f;

        //Make a 20MB buffer.
        samples = new byte[20000000];
        int sample_length = 0;

        //This block here simply loads the uncompressed data from the ogg file into a nice n' large 20MB buffer. If you want to use the same library as I've used, It's called DragonOgg (If you cant tell by the namespace)
        while (sample_length < samples.Length)
        {
            var bs = f.GetBufferSegment(4096); //Get ~4096 bytes (does not gurantee that 4096 bytes will be returned.
            if (bs.ReturnValue == 0)
                break; //End of file

            //Set the rate
            rate = bs.RateHz;

            //Display some loading info:
            Console.WriteLine("seconds: " + sample_length / rate);

            //It's stereo so we want half the data.
            int max = bs.ReturnValue / 2;

            //Buffer overflow care.
            if (samples.Length - sample_length < max)
                max = samples.Length - sample_length;

            //The copier.
            for (int j = 0; j < max; j++)
            {
                //I'm using j * 2 here because I know that the input audio is 8Bit Stereo, and we want just one mono channel. So we skip every second one.
                samples[sample_length + j] = bs.Buffer[j * 2];
            }

            sample_length += max;
            if (max == 0)
                break;
        }
        sampleCount = (sample_length - 1) / samplespacing + 1;
        frequencies = new System.Numerics.Complex[sampleCount][];
        for (int i = 0; i < sample_length; i += samplespacing)
        {
            Console.WriteLine("Sample---" + i + " / " + sample_length);
            System.Numerics.Complex[] sample;
            if (i + samplesize > sample_length)
                sample = new System.Numerics.Complex[sample_length - i];
            else
                sample = new System.Numerics.Complex[samplesize];
            for (int j = 0; j < sample.Length; j++)
            {
                sample[j] = (float)(samples[i + j] - 128) / 128.0f;
            }
            dft.BluesteinForward(sample, MathNet.Numerics.IntegralTransforms.FourierOptions.Default);
            frequencies[i / samplespacing] = sample;
        }

        //Perform the filters to the frequencies
        ExampleLowpass();

        //Make window kernel thingy
        float[] kernel = new float[samplesize / samplespacing * 2];
        for (int i=0; i<kernel.Length; i++)
        {
            kernel[i] = (float)((1-Math.Cos(2*Math.PI*i/(kernel.Length - 1)))/2);
        }

        //Apply window kernel thingy
        for (int i = 0; i < sample_length; i++)
        {
            int jstart = i / samplespacing - samplesize / samplespacing + 1;
            int jend = i / samplespacing;
            if (jstart < 0) jstart = 0;
            float ktotal = 0;
            float stotal = 0;
            for (int j = jstart; j <= jend; j++)
            {
                float kernelHere = 1.0f;
                if (jstart != jend)
                    kernelHere = kernel[(j - jstart) * kernel.Length / (jend + 1 - jstart)];
                int index = i - j * samplespacing;
                stotal += (float)frequencies[j][index].Real * kernelHere;
                ktotal += kernelHere;
            }
            if (ktotal != 0)
            {
                stotal /= ktotal;
                samples[i] = (byte)(stotal * 128 * 0.9f + 128);
            }
            else
            {
                Console.WriteLine("BAD " + jstart + " " + jend + " sec: " + ((float)i / rate));
                samples[i] = (byte)(stotal * 128 * 0.9f + 128);
            }
        }

        BinaryWriter bw = new BinaryWriter(File.OpenWrite("sound"));
        for (int i = 0; i < sample_length; i++)
        {
            bw.Write(samples[i]);
        }
        bw.Close();
    }
}

If you want to compile this, you'll need DragonOgg (http://sourceforge.net/projects/dragonogg/) and MathNet.Numerics (http://mathnetnumerics.codeplex.com/)

I hope it helps someone - I don't know about how StackOverflow licences by default, but this code as it stands is public domain.

On further thought, I decided I could achieve an approximated effect much easier by simply "blurring" the samples to get a basic low pass filter. A high pass filter can be made by subtracting the result of the low pass.