'impuls response code' von python nach matlab konvertieren

Question

Ich versuche einen Python Code in Matlab umzuwandeln. Die Code-Funktionen dienen dazu, eine Impulsantwort von dem ursprünglichen und dem aufgezeichneten Signal zu erhalten. Der Python-Code wird von Joseph erhalten. Ich glaube, der umgeschriebene Code ist derselbe wie der Python-Code. Die Impulsantwort von Python und Matlab ist jedoch unterschiedlich. Was habe ich hier falsch gemacht? Ich hoffe, ein paar Ratschläge und Kritik zu hören.

Vergleich der Impulsantwort von Python (links) und Matlab (rechts)

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 [sollte die Funktion wie Padarray im Python-Code]

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 [das gleiche wie filter20_20k (x, sr) funktionieren sollte]

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 

Python-Code von Joseph Ernest

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

Managed Joseph's impulseresponse.py python codes konvertieren in Matlab jetzt . Putting hier für zukünftige Möglichkeit von Leuten, die es verwenden.

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 und irrft.m

Gehen Sie zu diesem link. Der Dank geht an Dmitri Chubarov.