65 lines
1.6 KiB
Mathematica
65 lines
1.6 KiB
Mathematica
|
%
|
||
|
|
% getIR.m - extracts the IR from a logsweep output (input 20 to 20kHz, 10 s
|
||
|
|
% long, 2 s silence before and after)
|
||
|
|
%
|
||
|
|
% Copyright (C) Orastron srl unipersonale 2022
|
||
|
|
%
|
||
|
|
% This program is free software: you can redistribute it and/or modify
|
||
|
|
% it under the terms of the GNU General Public License as published by
|
||
|
|
% the Free Software Foundation, version 3 of the License.
|
||
|
|
%
|
||
|
|
% This program is distributed in the hope that it will be useful,
|
||
|
|
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||
|
|
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||
|
|
% GNU General Public License for more details.
|
||
|
|
%
|
||
|
|
% You should have received a copy of the GNU General Public License
|
||
|
|
% along with A-SID. If not, see <http://www.gnu.org/licenses/>.
|
||
|
|
%
|
||
|
|
% File author: Stefano D'Angelo
|
||
|
|
%
|
||
|
|
% Based on the theory in
|
||
|
|
%
|
||
|
|
% A.Farina, "Simultaneous measurement of impulse response and distortion with
|
||
|
|
% a swept-sine technique", 108th AES Convention, Paris, France, February 2000.
|
||
|
|
|
||
|
|
% y: logsweep output
|
||
|
|
% fs: sample rate (Hz)
|
||
|
|
% h: IR
|
||
|
|
function h = getIR(y, fs)
|
||
|
|
|
||
|
|
fstart=20;
|
||
|
|
fend=20e3;
|
||
|
|
silence=2;
|
||
|
|
T=10;
|
||
|
|
|
||
|
|
% sweep generation
|
||
|
|
t=[0:round(T*fs)-1];
|
||
|
|
t1=round(fs*silence);
|
||
|
|
t2=round(T*fs+silence*fs);
|
||
|
|
t3=round(T*fs+2*silence*fs);
|
||
|
|
w1 = 2*pi*fstart/fs;
|
||
|
|
w2 = 2*pi*fend/fs;
|
||
|
|
K = T*fs*w1/log(w2/w1);
|
||
|
|
L = T*fs/log(w2/w1);
|
||
|
|
sweep(1:t1-1)=0;
|
||
|
|
sweep(t1:t2-1)=sin(K.*(exp(t./L)-1));
|
||
|
|
sweep(t2:t3-1)=0;
|
||
|
|
|
||
|
|
% inverse filter
|
||
|
|
w = w1*exp(t./(T*fs).*log(w2/w1));
|
||
|
|
env = 10.^(-6/20.*log2(w./w1));
|
||
|
|
sweep_inv(1:t1-1)=0;
|
||
|
|
sweep_inv(t1:t2-1)=fliplr(sweep(t1:t2-1)).*env;
|
||
|
|
sweep_inv(t2:t3-1)=0;
|
||
|
|
|
||
|
|
% deconvolution
|
||
|
|
h=fftconv(sweep_inv,y);
|
||
|
|
|
||
|
|
% extract IR
|
||
|
|
[a,i]=max(abs(h));
|
||
|
|
i -= round(fs/2);
|
||
|
|
h=h(i:i+fs-1);
|
||
|
|
|
||
|
|
endfunction
|