Akustische Signalverarbeitung mit Python (2)

Dieser Artikel ist der 17. Tagesartikel von Python Advent Calendar. Ich werde über den Inhalt der akustischen Signalverarbeitung mit Python schreiben, den ich nicht in Artikel des letzten Jahres geschrieben habe.

Inhaltsverzeichnis

• Ich möchte die Abtastrate umrechnen
• Ich möchte eine 24-Bit-Audiodatei mit dem Wave-Modul lesen
• Ich möchte die Gruppenverzögerung mit scipy.signal.group_delay schätzen
• Ich möchte die Bildverarbeitung mit dem Schritttrick von numpy durchführen (ich möchte ein Spektrogramm zeichnen).
• Ich möchte ein TSP-Signal erzeugen
• Ich möchte Nullen und Pole visualisieren

Ich möchte die Abtastrate umwandeln

Um die Abtastrate umzuwandeln, müssen Sie, wenn die ursprüngliche Abtastrate und die Abtastrate, die Sie konvertieren möchten, ganzzahlige Verhältnisse sind, nur eine Aufwärts- / Abwärtsabtastung durchführen. Bei rationalen Verhältnissen verwenden Sie jedoch das minimale gemeinsame Vielfache jeder Abtastrate. Sie müssen auf die gewünschte Abtastrate hoch- und dann runtersampeln.

Hier betrachten wir die Konvertierung der Abtastrate von 48 kHz auf 44,1 kHz. Die Abtastratenumwandlung kann mit dem folgenden Code realisiert werden. Das Ziel sind die ersten 10 Sekunden einer 48-kHz-Tonquelle.

#!/usr/bin/env python
# vim:fileencoding=utf-8

from fractions import Fraction

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

import scipy.signal as sg

import soundfile as sf

if __name__ == '__main__':
    plt.close("all")

    fs_target = 44100
    cutoff_hz = 21000.0
    n_lpf = 4096

    sec = 10

    wav, fs_src = sf.read("../wav/souvenir_mono_16bit_48kHz.wav")
    wav_48kHz = wav[:fs_src * sec]

    frac = Fraction(fs_target, fs_src)  # 44100 / 48000

    up = frac.numerator  # 147
    down = frac.denominator  # 160

    # up sampling
    wav_up = np.zeros(np.alen(wav_48kHz) * up)
    wav_up[::up] = up * wav_48kHz
    fs_up = fs_src * up

    cutoff = cutoff_hz / (fs_up / 2.0)
    lpf = sg.firwin(n_lpf, cutoff)

    # filtering and down sampling
    wav_down = sg.lfilter(lpf, [1], wav_up)[n_lpf // 2::down]

    # write wave file
    sf.write("down.wav", wav_down, fs_target)

    # lowpass filter plot
    w, h = sg.freqz(lpf, a=1, worN=1024)
    f = fs_up * w / (2.0 * np.pi)
    fig = plt.figure(1)
    ax = fig.add_subplot(111)
    ax.semilogx(f, 20.0 * np.log10(np.abs(h)))
    ax.axvline(fs_target, color="r")
    ax.set_ylim([-80.0, 10.0])
    ax.set_xlim([3000.0, fs_target + 5000.0])
    ax.set_xlabel("frequency [Hz]")
    ax.set_ylabel("power [dB]")

    plt.show()

Da 44100/48000 = 147/160 ist, wird die ursprüngliche Tonquelle zuerst auf eine 147-fache Punktzahl hochgesampelt. Hier wird ein Puffer mit einer Größe von 147 Mal erstellt, und das ursprüngliche Signal wird in Intervallen von 147 Punkten in den Puffer eingesetzt. Infolgedessen werden nicht zugewiesene Punkte mit Nullen aufgefüllt. Entwerfen Sie dann den LPF so, dass er nur Frequenzen unterhalb der Nyquist-Frequenz von 44100 Hz enthält, und verwenden Sie lfilter, um den LPF anzuwenden, um das Band zu begrenzen. Der Punkt hier ist, dass die Nyquist-Frequenz während des Upsamplings 48000/2 = 24000 Hz beträgt, aber nach dem Downsampling 44100/2 = 22050 Hz zur Nyquist-Frequenz wird, sodass das während des Upsamplings angewendete Bandbegrenzungsfilter 22050 oder weniger beträgt. Es ist möglich, nur einen Filter zu verwenden, der wird. Nach der Bandbegrenzung wird nach Korrektur der Verzögerung aufgrund von LPF in Intervallen von 160 Punkten ein Downsampling durchgeführt. Als Ergebnis wurde eine Abtastratenumwandlung realisiert.

Die tatsächliche Audio-Wellenform und das Spektrogramm sind wie folgt. Die oberen beiden befinden sich nach der 44,1-kHz-Konvertierung und die unteren beiden befinden sich vor der Konvertierung.

Ich möchte eine 24-Bit-Audiodatei mit dem Wave-Modul ## lesen

Wenn Sie eine 16-Bit-Audiodatei mit einem Standard-Wave-Modul lesen möchten, können Sie sie mit wave.open lesen und mit readframes so viel lesen, wie Sie benötigen, und sie mit np.frombuffer in int16 konvertieren. Bei einer 24-Bit-Soundquelle können Sie jedoch kein 24-Bit mit frombuffer angeben. Daher müssen Sie es selbst lesen. Hier wird, wie im folgenden Code gezeigt, beim Lesen von jeweils 3 Bytes mit dem Entpacken des Strukturmoduls die 24-Bit-Tonquelle durch Packen von 0 und Entpacken als int32 gelesen.

Die zu lesende Tonquelle ist übrigens e-onkyo sample sound source.

#!/usr/bin/env python
# vim:fileencoding=utf-8

import numpy as np
import matplotlib.pyplot as plt
import wave

from struct import unpack

if __name__ == '__main__':
    plt.close("all")

    fname = "../wav/souvenir.wav"
    fp = wave.open(fname, "r")

    nframe = fp.getnframes()
    nchan = fp.getnchannels()
    nbyte = fp.getsampwidth()
    fs = fp.getframerate()

    print("frame:{0}, "
          "channel:{1}, "
          "bytewidth:{2}, "
          "fs:{3}".format(nframe, nchan, nbyte, fs))

    buf = fp.readframes(nframe * nchan)
    fp.close()

    read_sec = 40
    read_sample = read_sec * nchan * fs
    print("read {0} second (= {1} frame)...".format(read_sec,
                                                    read_sample))

    #Durch Packen von 0 in das niedrigste Bit und Auspacken von int
    #Extrahieren Sie den Wert mit dem 24-Bit-Wert als 32-Bit-Int
    #  (<Ich nehme einen kleinen andian int Wert an)
    #Auspacken gibt Tupel zurück[0]ich nehme
    unpacked_buf = [unpack("<i",
                           bytearray([0]) + buf[nbyte * i:nbyte * (i + 1)])[0]
                    for i in range(read_sample)]

    #ndarray
    ndarr_buf = np.array(unpacked_buf)

    # -1.0〜1.Normalisiert auf 0
    float_buf = np.where(ndarr_buf > 0,
                         ndarr_buf / (2.0 ** 31 - 1),
                         ndarr_buf / (2.0 ** 31))

    #Interleave lösen(Für Stereo-Tonquelle)
    wav_l = float_buf[::2]
    wav_r = float_buf[1::2]
    time = np.arange(np.alen(wav_l)) / fs

    # plot
    fig = plt.figure(1)

    ax = fig.add_subplot(2, 1, 1)
    ax.plot(time, wav_l)
    ax.set_xlabel("time [pt]")
    ax.set_ylabel("amplitude")
    ax.set_title("left channel")
    ax.grid()

    ax = fig.add_subplot(2, 1, 2)
    ax.plot(time, wav_r)
    ax.set_xlabel("time [pt]")
    ax.set_ylabel("amplitude")
    ax.set_title("right channel")
    ax.grid()

    plt.show()

Ich möchte die Gruppenverzögerung mit scipy.signal.groupdelay ## schätzen

[Dieser Artikel] aus dem letztjährigen Artikel (http://qiita.com/wrist/items/5759f894303e4364ebfd#%E7%BE%A4%E9%81%85%E5%BB%B6%E3%82%92%E8% Ich habe es nicht bemerkt, als ich A8% 88% E7% AE% 97% E3% 81% 97% E3% 81% 9F% E3% 81% 84 geschrieben habe), aber ich habe es nicht bemerkt, aber die [group_delay-Methode] zum Finden der Gruppenverzögerung ( https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.signal.group_delay.html#scipy.signal.group_delay) wurde seit Version 0.16.0 von scipy.signal hinzugefügt. Es war. Dies kann verwendet werden, um die Gruppenverzögerung eines digitalen Filters zu bestimmen.

#!/usr/bin/env python
# vim:fileencoding=utf-8

import numpy as np
import matplotlib.pyplot as plt

import scipy.signal as sg


def allpass_filter(freq, r=0.9, fs=16000.0):
    omega = 2.0 * np.pi * freq / fs

    b = [r ** 2, - 2.0 * r * np.cos(omega), 1.0]
    a = [1.0, -2.0 * r * np.cos(omega), r ** 2]

    return (b, a)

if __name__ == '__main__':
    plt.close("all")

    fs = 16000.0
    N = 1024

    b, a = allpass_filter(1000.0)
    w, h = sg.freqz(b, a, worN=N)

    f = w * fs / (2.0 * np.pi)
    _, gd = sg.group_delay((b, a), w=w)

    fig = plt.figure(1)
    ax = fig.add_subplot(3, 1, 1)
    ax.semilogx(f, np.abs(h))
    ax.grid()
    ax.set_xlabel("frequency [Hz]")
    ax.set_ylabel("amplitude")
    ax.set_xlim([10.0, fs / 2.0])

    ax = fig.add_subplot(3, 1, 2)
    ax.semilogx(f, np.angle(h))
    ax.grid()
    ax.set_xlabel("frequency [Hz]")
    ax.set_ylabel("phase [rad]")
    ax.set_xlim([10.0, fs / 2.0])

    ax = fig.add_subplot(3, 1, 3)
    ax.semilogx(f, gd)
    ax.grid()
    ax.set_xlabel("frequency [Hz]")
    ax.set_ylabel("group delay [pt]")
    ax.set_xlim([10.0, fs / 2.0])

    plt.show()

Ich möchte die Frame-Verarbeitung mit dem Schritttrick von numpy durchführen (ich möchte ein Spektrogramm zeichnen)

Unter Verwendung des in den Rezepten 4.6 und 4.7 des IPython Data Science Cookbook eingeführten Schritttricks wird ein Puffer in kurze Frames unterteilt. Sie können eine Frame-Verarbeitung effizient durchführen, auf die wiederholt zugegriffen wird, ohne eine Kopie des Frames zu generieren. Mit as_strided unter numpy.lib.stride_tricks können Sie die Methode zur Berechnung des Spektrogramms wie folgt schreiben.

#!/usr/bin/env python
#vim:fileencoding=utf-8

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import matplotlib.cm as cm

import scipy.fftpack as fft
from numpy.lib.stride_tricks import as_strided

import soundfile as sf

def spectrogram(xs, nfft=1024, noverlap=2, win=np.hanning):

    nsample = np.alen(xs)
    nframe = int(nsample * (noverlap / float(nfft))) - noverlap
    nbyte = xs.dtype.itemsize
    shift = (nfft // noverlap)

    xs_split = as_strided(xs, (nframe, nfft), (nbyte * shift, nbyte))
    window = win(nfft)

    Xs = []
    print(nsample, nframe, nfft, shift, nbyte, nfft * nframe)
    for frame in xs_split:
        Xs.append(20.0 * np.log10(np.abs(fft.fft(window * frame, nfft))[:nfft//2 + 1]))

    return np.array(Xs)


if __name__ == '__main__':
    plt.close("all")

    wav, fs = sf.read("../wav/souvenir_mono_16bit_48kHz.wav")

    sec = 20
    nread = fs * sec

    nfft = 1024
    spec = spectrogram(wav[:nread], nfft=nfft, noverlap=2)

    # plot
    fig = plt.figure(1)
    ax = fig.add_subplot(211)
    ax.plot(wav[:nread])
    ax.grid()

    ax = fig.add_subplot(212)
    ax.pcolor(np.arange(spec.shape[0]),
              np.arange(nfft//2 + 1) * (fs / nfft),
              spec.T,
              cmap=cm.jet)
    ax.set_xlim([0, spec.shape[0]])
    ax.set_ylim([50, fs//2])
    ax.set(yscale="log")

    plt.show()

Schritt ist ein Wert, der das Byte-Intervall beim Zugriff auf einen zusammenhängenden Bereich darstellt. Wenn Sie dies ändern, können Sie steuern, wie viele Bytes Sie ausführen, wenn Sie einen Zeilen- oder Spaltenindex vorrücken.

as_strided ist eine Methode, um auf den Puffer zuzugreifen, wobei dieser Schritt pseudo geändert wird. Das erste Argument ist der Zugriffszielpuffer, das zweite Argument ist die Form des Puffers, der durch Ändern des Schritts erstellt wird, und das dritte Argument ist der Schritt, den Sie ändern möchten. Wenn Sie nach dem Ändern des Puffers des ersten Arguments auf den Schritt des dritten Arguments in Pseudo-Weise zugreifen, wird die Operation so ausgeführt, dass ein Alias zurückgegeben wird, auf den als Puffer zugegriffen werden kann, der durch die Form des zweiten Arguments angegeben wird.

Das Format dieses Schrittes ist (Byte-Abstand beim Vorrücken von Zeilen, Byte-Abstand beim Vorrücken von Spalten). Im Fall des obigen Skripts ist das ursprüngliche Argument xs ein eindimensionales Array und sein Schritt ist (8,). Wenn Sie also den Index um 1 vorrücken, wird er um 8 Bytes vorgerückt (dtype ist float64). Wenn dieser Schritt beispielsweise durch as_stride in (512 * 8, 8) geändert wird, kann xs, das ursprünglich ein eindimensionales Array war, als zweidimensionales Array xs_split betrachtet werden. Wenn Sie wiederholt mit einer for-Schleife auf xs_split zugreifen, werden Sie jedes Mal, wenn Sie darauf zugreifen, in die Zeilenrichtung von xs_split vorrücken, dh Sie werden 512 * 8 Bytes xs vorrücken, sodass Sie darauf zugreifen können, indem Sie xs um 512pt verschieben. Auf diese Weise kann die Rahmenverarbeitung, die andernfalls durch Ausschneiden eines eindimensionalen Arrays in mehrere Kurzzeitrahmen verarbeitet würde, jetzt durch Schritttricks ausgeführt werden, ohne zum Zeitpunkt des Ausschneidens eine Kopie zu erzeugen. Ich werde.

Beachten Sie, dass dieser Schritttrick den Zugriff außerhalb des Bereichs des internen Puffers erleichtert. Wenn Sie also einen Fehler machen, greifen Sie häufig auf den undefinierten Bereich zu und fallen.

Ich möchte ein TSP-Signal ## erzeugen

Wenn ich einen Sweep-Impuls (TSP-Signal) entwerfe, der für die akustische Messung verwendet wird, denke ich, dass er häufig durch IFFTing erhalten wird, was im Frequenzbereich entworfen wurde, und in den Zeitbereich zurückgeführt wird. Es kann mit dem folgenden Code generiert werden. N_exp und m_exp, die zur effektiven Länge beitragen, müssen angepasst werden.

#!/usr/bin/env python
# vim:fileencoding=utf-8

import numpy as np
import matplotlib.pyplot as plt

import scipy.fftpack as fft

import soundfile as sf

if __name__ == '__main__':
    plt.close("all")

    # parameters
    N_exp = 16
    m_exp = 2
    nrepeat = 5
    fs = 48000
    gain = 100.0

    N = 2 ** N_exp
    m = N // (2 ** m_exp)  # (J=2m)
    a = ((m * np.pi) * (2.0 / N) ** 2)

    tsp_freqs = np.zeros(N, dtype=np.complex128)
    tsp_freqs[:(N // 2) + 1] = np.exp(-1j * a * (np.arange((N // 2) + 1) ** 2))
    tsp_freqs[(N // 2) + 1:] = np.conj(tsp_freqs[1:(N // 2)][::-1])

    # ifft and real
    tsp = np.real(fft.ifft(tsp_freqs, N))

    # roll
    tsp = gain * np.roll(tsp, (N // 2) - m)

    # repeat
    tsp_repeat = np.r_[np.tile(tsp, nrepeat), np.zeros(N)]

    # write
    sf.write("tsp.wav", tsp_repeat, fs)

    fig = plt.figure(1)
    ax = fig.add_subplot(211)
    ax.plot(tsp)
    ax = fig.add_subplot(212)
    ax.plot(tsp_repeat)

    plt.show()

Die dadurch generierte WAV-Datei wird wie unten gezeigt als Audacity angezeigt.

Das TSP-Signal wird in den folgenden grundlegenden Schulungsmaterialien für die Messung der Impulsantwort von Professor Kaneda ausführlich beschrieben. http://www.asp.c.dendai.ac.jp/ASP/IRseminor2016.pdf

Ich möchte die Nullen und Pole visualisieren

Sie können den Zeitfilterkoeffizienten mit tf2zpk in Null und Pol konvertieren. Um dies zu veranschaulichen, können Sie mit add_patch wie unten gezeigt einen Kreis zeichnen und dann ohne Linien und mit Markierungen zeichnen. .. Ich habe es nach dem Einfügen der Figur bemerkt, aber da es sich um einen FIR-Filter handelt, gibt es keinen Pol, aber bitte ignorieren Sie ihn.

#!/usr/bin/env python
#vim:fileencoding=utf-8

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

import scipy.signal as sg

if __name__ == '__main__':
    plt.close("all")

    b = sg.firwin(11, 0.2, window="han")
    z, p, k = sg.tf2zpk(b, [1])

    fig = plt.figure(1, figsize=(8, 8))
    ax = fig.add_subplot(111)

    ax.add_patch(plt.Circle((0.0, 0.0), 1.0, fc="white"))
    ax.plot(np.real(z), np.imag(z), "o", mfc="white")
    ax.plot(np.real(p), np.imag(p), "x", mfc="white")

    ax.grid()
    ax.set_xlim([-1.5, 1.5])
    ax.set_ylim([-1.5, 1.5])

    plt.show()

Andere

Ich wollte auch den folgenden Inhalt schreiben, hatte aber nicht genug Zeit, daher möchte ich ihn erneut hinzufügen oder in einen anderen Artikel schreiben.

• Ich möchte die Impulsantwort aus den Eingangs- / Ausgangssignalen abschätzen
• Konvertierung des Digitalfilters in das SOS-Format
• Ich möchte die Signalhüllkurve mithilfe der Hilbert-Transformation ermitteln

In Bezug auf die Hilbert-Konvertierung habe ich nur die folgende Erklärung geschrieben, daher werde ich sie vorerst veröffentlichen.

Um die Hüllkurve zu finden, die die maximale Amplitude des Signals reibungslos verfolgt, gibt es ein Verfahren zum Multiplizieren des Signals mit dem Absolutwert oder Quadratwert mit LPF, aber auch das Verfahren zum Ermitteln aus dem durch Hilbert-Umwandlung erhaltenen Analysesignal. Ist gewesen

Referenz https://jp.mathworks.com/help/dsp/examples/envelope-detection.html