将“获取脉冲响应代码”从 python 转换为 matlab答案

【问题标题】：Convert 'obtain impulse response code' from python to matlab将“获取脉冲响应代码”从 python 转换为 matlab
【发布时间】：2017-08-20 22:18:04
【问题描述】：

我正在尝试将 python 代码转换为 Matlab。代码功能是从原始信号和记录信号中获得脉冲响应。 python代码来自Joseph。我相信重写的代码和python代码是一样的。但是从 python 和 Matlab 获得的脉冲响应是不同的。我在这里做错了什么？希望听到一些建议和批评。

python(左)和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 [应该与 filter20_20k(x,sr) 功能相同]

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

Joseph Ernest 编写的 Python 代码

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)

【问题讨论】：

在将代码转换为另一种语言时发现错误的一个好方法是评估中间结果或独立测试部分代码
@m7913d 根据您的建议，我尝试通过打印值来查看值。但由于阵列规模大（960 万），它被卡住了。所以我尝试将它绘制在图表上，但出现此错误OverflowError: In draw_path: Exceeded cell block limit。我在 matlab 但不是 python 上查看结果没有问题。我应该如何改变我的尝试？
用一个小例子测试你的代码。

标签： python matlab

【解决方案1】：

现在设法将Joseph's impulseresponse.py python codes 转换为matlab。放在这里是为了人们将来使用它的可能性。

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。

【讨论】：