Apply frequency sampling filter on audio file

Question

I am trying to apply a double-bandstop filter using frequency sampling from a given data file. The method I am using is as follows

1. Inverse fft (ifft) of the given frequency sampling file

2. Circular shifting of the real values given from step 1

3. Perform an fft of the result from step 2

4. Apply the filter on the audio file using convolution. (frequency domain)

The problem is that the bandstop frequencies (925Hz & 2090Hz) still exist. Is there any problem with my code or do I miss something?

[wave,fs]=audioread('audio.WAV'); 
data=importdata ('freqSampling.txt')

y=(ifft(data,401))
x=real (y)

r=circshift (x,200)

f=fft (r,4096);

new_sound=conv (wave, f)
sound(new_sound,fs,16);

Can anyone help me with that?

Answer 1

The first two steps

1. Inverse fft (ifft) of the given frequency sampling file
2. Circular shifting of the real values given from step 1

should correctly give you the time-domain coefficients of a filter constructed from your specifications, provided the freqSampling.txt file correctly specifies the desired full double-sided spectrum (see "Validating specifications" below). The number of ifft points may also need to be adjusted/increased if the frequency specification contains steep transients. Performing the convolution in the frequency domain as you indicated in step 4 however does not correspond to a typical filtering operation, but would rather be equivalent to a time-domain multiplication of the two signals.

Time-domain filtering

From the result of step 2 you could filter your wave data directly in the time-domain using either conv:

new_sound = conv(r, wave);

or filter:

new_sound = filter(r, 1, wave);

conv would be giving you the full length(wave)+length(r)-1 convolution, whereas filter is a more signal processing oriented function returning the first length(wave) samples of the convolution (and can also handle recursive filters which conv does not support directly).

Frequency-domain filtering

Alternatively to perform the filtering in the frequency-domain you would

3. Perform an FFT of the result from step 2, using a size that is at least length(r)+length(wave)-1
4. Perform an FFT of the wave data, using the same size as in step 3
5. Apply the filter on the audio file using multiplication (frequency domain).
6. Compute the inverse FFT (ifft) of the result from step 5
7. Take the real part of the result from step 6 (technically not required, but often needed due to small numerical round-off errors)

This can be implemented using the following:

N = length(wave)+length(r)-1;
wave_fd = fft(wave, N);              % step 3
filter_fd = fft(r, N);               % step 4

filtered_fd = wave_fd .* filter_fd;  % step 5
new_sound = real(ifft(filtered_fd)); % step 6 & 7

Note that you could also perform this frequency-domain filtering operation in smaller chunks using the overlap-add method.

Validating specifications

Based on your comments, the data imported from the freqSampling.txt file could be reconstructed with:

N = 401;
data = ones(N,1);
data(19:23) = [2 1 0 1 2]/3;
data(51:56) = [2 1 0 1 2]/3;
data(N-[2:(N+1)/2]+2) = data([2:(N+1)/2]);

To validate that this would filter the desired frequencies, we can plot this specification as a function of frequency. To do this we would need the sampling rate used (fs), which seem to be 22050 according to your graph. You should then be able to plot these with:

hold off; plot([0:N-1]*fs/N, data);
hold on;  plot([925 925;2090 2090]', [0 1.2;0 1.2]', 'k:');
axis([0 3000 0 1.2]);
xlabel('Frequency (Hz)');
ylabel('Amplitude');
legend('Specs', 'Tones');

Which should give a plot that looks like:

Based on this, it seems that the specifications do not provide any attenuation at the tone frequencies. A better fit could be constructed with:

N = 401;
data = ones(N,1);
data(round(925*N/fs)+1+[-2:2]) = [2 1 0 1 2]/3;  % data([16:20])
data(round(2090*N/fs)+1+[-2:2]) = [2 1 0 1 2]/3; % data([37:41])
data(N-[2:(N+1)/2]+2) = data([2:(N+1)/2]);

Yielding a validation plot that looks like:

P.S.: based on your signal's frequency-domain graph, it would seem that the second tone is closer to 2600Hz rather than the indicated 2090Hz.