Pythonからmatlabに 'impulse response codeを取得する'

Question

PythonコードをMatlabに変換しようとしています。コード機能は、元の信号と記録された信号からインパルス応答を得ることです。 PythonコードはJosephから得られます。私は、書き直されたコードは、Pythonコードと同じだと考えています。しかし、PythonとMatlabから得られるインパルス応答は異なります。私はここで何をやったのですか？いくつかのアドバイスや批判を聞いて欲しい。パイソン（左）および（右）MATLAB

impulseresponse.mからインパルス応答の

比較[Matlabの]

[a,fs] = audioread('sweep.wav'); % sweep 
[b,fs] = audioread('rec.wav'); % rec 

a = pad(a,fs*50,fs*10); 
b = pad(b,fs*50,fs*10); 
[m,n] = size(b); 
h = zeros(m,n); 

for chan = 1:2 
    b1 = b(:,chan); 
    b1 = filter20_20k(b1,fs); 
    ffta = fft(a); 
    fftb = fft(b1); 
    ffth = fftb ./ ffta; 
    h1 = ifft(ffth); 
    h1 = filter20_20k(h1,fs); 
    h(:,chan) = h1; 
end 

h = h(1:10*fs,:); 
dB = 40; 
h = h*power(10,dB*1.0/20); 
h(:,1) = h(:,1) ./ max(abs(h(:,1))); % normalizing impulse response 
h(:,2) = h(:,2) ./ max(abs(h(:,2))); 
figure; 
plot(h) 

pad.m [機能すべきPythonコードのpadarrayと同じ]

function y = pad(data, t_full, t_pre) 
[row_dim,col_dim] = size(data); 
t_post = t_full - row_dim - t_pre; 
if t_post > 0 
    if col_dim == 1 
     y = [zeros(t_pre,1);data;zeros(t_post,1)]; 
    else 
     y1 = [zeros(t_pre,1);data(:,1);zeros(t_post,1)]; 
     y2 = [zeros(t_pre,1);data(:,2);zeros(t_post,1)]; 
     y = [y1,y2]; 
    end 
else 
    if col_dim == 1 
     y = [zeros(t_pre,1);data(1:t_full - t_pre,1)]; 
    else 
     y1 = [zeros(t_pre,1);data(1:t_full - t_pre,1)]; 
     y2 = [zeros(t_pre,1);data(1:t_full - t_pre,2)]; 
     y = [y1,y2]; 
    end 
end 
end 

filter20_20k.mジョセフアーネスト今てMATLABにJoseph's impulseresponse.py python codesを変換するために管理

SWEEPFILE = 'sweep.wav' 
RECFILE = 'rec.wav' 
OUTPUTFILE = 'IR.wav' 

import wavfile 
from scipy import signal 
import numpy as np 

def ratio(dB): 
    return np.power(10, dB * 1.0/20) 

def padarray(A, length, before=0): 
    t = length - len(A) - before 
    if t > 0: 
     width = (before, t) if A.ndim == 1 else ([before, t], [0, 0]) 
     return np.pad(A, pad_width=width, mode='constant') 
    else: 
     width = (before, 0) if A.ndim == 1 else ([before, 0], [0, 0]) 
     return np.pad(A[:length - before], pad_width=width, mode='constant') 

def filter20_20k(x, sr): # filters everything outside out 20 - 20000 Hz 
    nyq = 0.5 * sr 
    sos = signal.butter(5, [20.0/nyq, 20000.0/nyq], btype='band', output='sos') 
    return signal.sosfilt(sos, x) 

sr, a, br = wavfile.read(SWEEPFILE, normalized=True) 
sr, b, br = wavfile.read(RECFILE, normalized=True) 

a = padarray(a, sr*50, before=sr*10) 
b = padarray(b, sr*50, before=sr*10) 
h = np.zeros_like(b) 

for chan in [0, 1]: 
    b1 = b[:,chan] 

    b1 = filter20_20k(b1, sr) 
    ffta = np.fft.rfft(a) 
    fftb = np.fft.rfft(b1) 
    ffth = fftb/ffta 
    h1 = np.fft.irfft(ffth) 
    h1 = filter20_20k(h1, sr) 

    h[:,chan] = h1 

h = h[:10 * sr,:] 
h *= ratio(dB=40) 

wavfile.write(OUTPUTFILE, sr, h, normalized=True, bitrate=24) 

Answer 1

によって

function y = filter20_20k(x,fs) 
nyq = 0.5*fs; 
[z,p,k] = butter(5,[20.0/nyq,20000.0/nyq],'bandpass'); 
sos = zp2sos(z,p,k); 
y = sosfilt(sos,x); 
end 

[（X、SR）filter20_20kと同様に機能するはず] Pythonコード。それを使用する人々の将来の可能性のためにここに置く。

impulseresponse.m

[a,fs] = audioread('sweep.wav'); % sweep 
[b,fs] = audioread('rec.wav'); % rec 

a = pad(a,fs*50,fs*10); 
b = pad(b,fs*50,fs*10); 
[m,n] = size(b); 
h = zeros(m,n); 

for chan = 1:2 
    b1 = b(:,chan); 
    b1 = filter20_20k(b1,fs); 
    ffta = rfft(a); 
    fftb = rfft(b1); 
    ffth = fftb ./ ffta; 
    h1 = irfft(ffth); 
    h1 = filter20_20k(h1,fs); 
    h(:,chan) = h1; 
end 

h = h(1:10*fs,:); 
dB = 40; 
h = h*power(10,dB*1.0/20); 

audiowrite('heck.wav',h,fs,'BitsPerSample',24); 
[y,fs] = audioread('heck.wav'); 
figure; 
plot(y) 

pad.m

function y = pad(data, t_full, t_pre) 
[row_dim,col_dim] = size(data); 
t_post = t_full - row_dim - t_pre; 
if t_post > 0 
    if col_dim == 1 
     y = [zeros(t_pre,1);data;zeros(t_post,1)]; 
    else 
     y1 = [zeros(t_pre,1);data(:,1);zeros(t_post,1)]; 
     y2 = [zeros(t_pre,1);data(:,2);zeros(t_post,1)]; 
     y = [y1,y2]; 
    end 
else 
    if col_dim == 1 
     y = [zeros(t_pre,1);data(1:t_full - t_pre,1)]; 
    else 
     y1 = [zeros(t_pre,1);data(1:t_full - t_pre,1)]; 
     y2 = [zeros(t_pre,1);data(1:t_full - t_pre,1)]; 
     y = [y1,y2]; 
    end 
end 
end 

filter20_20k.m

function y = filter20_20k(x,fs) 
nyq = 0.5*fs; 
[z,p,k] = butter(5,[20.0/nyq,20000.0/nyq],'bandpass'); 
sos = zp2sos(z,p,k); 
y = sosfilt(sos,x); 
end 

rfft.mとirrft.m

このlinkに移動します。クレジットはDmitri Chubarovに払い戻されます。