sr = 16000;

dt = 1/sr;

n_fft = 1024;

df = sr/n_fft;

Firstly, we plot the source signal and impulse response in both time domain and frequency domain.

x = load('src.txt');

plot(x(1:100));

title('time domain source signal');

xticklabels((0:10:100)*dt*1000);

xlabel('(ms)');

spec_x = fft(x,n_fft);

plot(abs(spec_x(1:n_fft/2+1)));

title('spectrum of source signal');

xlim([0 n_fft/2+1]);

xticklabels(0:sr/2/10:sr/2);

y = load('rec1.txt');

plot(y);

title('time domain impulse response');

xticklabels((0:length(y)/10:length(y))*dt*1000);

xlabel('(ms)');

spec_y = fft(y,n_fft);

plot(abs(spec_y(1:n_fft/2+1)));

title('spectrum of impulse response');

xlim([0 n_fft/2+1]);

xticklabels(0:sr/2/10:sr/2);

When we observe the transfer function we may find that the high frequency part (above 4000 Hz) of the transfer function is invalid, because the source signal even not have this frequency component. We need design a low pass filter.

tf = spec_y./spec_x;

plot(abs(tf(1:n_fft/2+1)));

title('transfer function');

xlim([0 n_fft/2+1]);

xticklabels(0:sr/2/10:sr/2);

fcut = 4000;

b = fir1(32,fcut/(sr/2));

tf_lowpass = tf.*fft(b,n_fft)';

plot(abs(tf_lowpass(1:n_fft/2+1)));

title('bandpass transfer function')

xlim([0 n_fft/2+1]);

xticklabels(0:sr/2/10:sr/2);

We use `fir1` to design a fir filter. For a low pass filter, we need to set the cut frequency, and pass the relative cut frequency as the second parameter to the fir function. The first parameter is the length of filter, the output `b` is a time domain filter.

ir = ifft(tf_lowpass);

plot(ir)

s = audioread('fbcut_16k.wav');

s = s./max(s)/2;

s_reverb = fftfilt(ir,s);

s_reverb = s_reverb./max(s_reverb)/2;

subplot(211);

plot(s);

title('original sound');

subplot(212);

plot(s_reverb);

title('convovled sound');