How to modify and then inverse complex data that went through FFT algorithm?

Question

I wrote a low-pass filter based around the FFT algorithm.

First I process the input data using forward FFT, then I decrease the "volume" of parts of the resulting spectrum, then I put the data into the inverse FFT and finally I normalize the data like this:

fn lowpass_filter(data: &[f32], sampling_rate: f32, cutoff_frequency: f32) -> Vec<f32> {
    let len = data.len();

    // step 1:
    let mut spectrum = fft::forward(&data);

    // step 2:
    let start: usize = (len as f32 * (cutoff_frequency / sampling_rate * 2.0)) as usize;
    let diff = len - start;

    // what to add to multiplikator in each iteration
    let step: f32 = PI / diff as f32;

    // reduce volume of frequencies after cutoff_frequency in spectrum
    let mut multiplikator: f32 = 0.0;
    for i in start..len {
        let mul = (multiplikator.cos() + 1.0) / 2.0;
        spectrum[i] *= mul;
        multiplikator += step;
     }

    // step 3:
    let data = fft::inverse(&spectrum);

    // step 4:
    fft::normalize(&data, true)
}

When only doing step 1, 3 and 4 it works but the only problem there is that after normalization of the inversed data, I only get the absolute value of the data back so a 65Hz sinus wave looks like this:

The main problem I am facing tho, is that I do not know how to reduce the volume of specific frequencies in the spectrum.

When reducing said frequencies like in step 3 the visualization of a 65Hz sinus wave put through that lowpass filter with the cutoff_frequency set to 100.0Hz looks like this:

What did I do wrong here?

More info about the defined functions:

use rustfft::FftPlanner;
pub use rustfft::num_complex::Complex;

pub fn forward(data: &[f32]) -> Vec<Complex<f32>> {
    let length = data.len();

    // conversion to complex numbers
    let mut buffer: Vec<Complex<f32>> = Vec::new();
    for d in data {
        buffer.push(Complex{re: *d, im: 0.0});
    }

    // creates a planner
    let mut planner = FftPlanner::<f32>::new();

    // creates a FFT
    let fft = planner.plan_fft_forward(length);

    //input.append(&mut data.to_vec());

    fft.process(&mut buffer);

    buffer
}

pub fn inverse(data: &[Complex<f32>]) -> Vec<Complex<f32>> {
    let length = data.len();

    let mut  data = data.to_vec();

    // creates a planner
    let mut planner = FftPlanner::<f32>::new();

    // creates a FFT
    let fft = planner.plan_fft_inverse(length);


    fft.process(&mut data);

    data.to_vec()
}

pub fn normalize(data: &[Complex<f32>], normalize: bool) -> Vec<f32> {
    let len: f32 = data.len() as f32;
    let norm = data
        .iter()
        .map(|x| x.norm() / len)
        .collect();

    norm
}

I am using rustfft for the actual processing.

Answer 1

There are two issues with your code:

1. Since you want to process real data (i.e. data whose imaginary part is 0), the output of the forward FFT is symmetrical and you need to apply the same coefficient to matching frequencies (i.e. spectrum[i] and spectrum[spectrum.len() - i], remembering that spectrum[0] stands alone, and so does spectrum[spectrum.len()/2] if the length is even).

2. If the frequencies are symmetrical, the result of the inverse FFT should be real (i.e. the imaginary part should be 0 ± small rounding errors). Therefore your normalize function should use x.re instead of x.norm(), which will allow it to retain its sign.

Answer 2

After fixing these issues and adding a cutoff_end_freq to the lowpass_filter function the code looks like this:

pub fn lowpass_filter(data: &[f32], sampling_rate: f32, cutoff_start_freq: f32, cutoff_end_freq: f32) -> Vec<f32> {
    let len = data.len();
    let spectrum_len = len / 2;

    let mut spectrum = fft::forward(&data);
    assert!(len == spectrum.len());

    let start: usize = (spectrum_len as f32 * (cutoff_start_freq / sampling_rate * 2.0)) as usize;
    let end: usize = (spectrum_len as f32 * (cutoff_end_freq / sampling_rate * 2.0)) as usize;
    let diff = end - start;

    // what to add to multiplikator in each iteration
    let step: f32 = PI / diff as f32;

    let mut multiplikator: f32 = 0.0;
    for i in start..=end {
        let mul = (multiplikator.cos() + 1.0) / 2.0;
        spectrum[i] *= mul;
        spectrum[len - i - 1] *= mul;

        multiplikator += step;
    }
    for i in end..spectrum_len {
        spectrum[i] *= 0.0;
        spectrum[len - i - 1] *= 0.0;
    }

    let data = fft::inverse(&spectrum);

    fft::get_real(&data)
}

The normalize function now renamed to get_real:

pub fn get_real(data: &[Complex<f32>]) -> Vec<f32> {
    let len: f32 = data.len() as f32;
    let norm = data
        .iter()
        .map(|x| x.re / len)
        .collect();

    norm
}

Now it works very well, but there are still alternating "spikes" on each end when visualizing a frequency outside of the threshold:

I simply put a Hanning windowing effect over it because I heard that the FFT expects the data to be continuous. And it works, with the drawback that a lot of information on the sides get lost.

Are there better methods to circumvent this problem or is there still something wrong with the lowpass?