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?
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.
There are two issues with your code:
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).
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.