-
Notifications
You must be signed in to change notification settings - Fork 12
/
convolutionReverb.m
80 lines (64 loc) · 1.93 KB
/
convolutionReverb.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
%% Load the input file and impulse response.
inputFilename = 'audio/singing.mp3';
% inputFilename = 'mozart.mp3';
irFilename = 'audio/stalbans_a_mono.wav';
[signal, fs] = audioread(inputFilename);
[impulseResponse, irFs] = audioread(irFilename);
assert(fs == irFs);
% Truncate for faster processing.
impulseResponse = impulseResponse(1:50000);
figure;
plot(impulseResponse);
%% Apply reverb using convolution.
reverberated = conv(signal, impulseResponse);
%% Plot and play the original and reverberated sound.
paddedSignal = [signal; zeros(length(impulseResponse) - 1, 1)];
%subplot(2,1,1);
figure;
hold on;
plot(reverberated);
plot(paddedSignal);
hold off;
legend('Reverberated', 'Original (dry)');
soundsc([signal; reverberated], fs);
%% Inverse filtering - frequency domain
nfft = length(reverberated);
hf = fft(impulseResponse, nfft);
spectrum = fft(reverberated);
spectrum = spectrum ./ hf;
% To reduce noise amplification:
% spectrum = spectrum .* conj(hf) ./ (abs(hf).^2 + 1e-2);
% See: https://blogs.mathworks.com/steve/2007/08/13/image-deblurring-introduction/
inverseFiltered = real(ifft(spectrum));
soundsc(inverseFiltered, fs);
%% Plot the signals frequency responses
figure;
freqz(reverberated);
title('Reverberated');
figure;
freqz(inverseFiltered);
title('Inverse filtered');
figure;
freqz(impulseResponse);
title('Impulse response');
%% Inverse filtering - time domain
% Doesn't work.
[dry, r] = deconv(reverberated, impulseResponse);
figure;
plot(dry);
soundsc(dry, fs);
% Why does this not work?
% In the beginning, it resembles an periodic signal,
% but explodes very fast.
figure;
freqz(impulseResponse);
% figure;
% zplane(impulseResponse);
inverseFiltered = filter(1, impulseResponse, reverberated);
freqz(1, impulseResponse);
% The inverse filter is not stable because the original filter
% is not minimal phase. The original filter has zeros on or outside of the
% unit circle.
isstable(1, impulseResponse)
figure;
plot(inverseFiltered);