Traitement du signal acoustique avec Python (2)

Cet article est l'article du 17ème jour du Calendrier de l'Avent Python. J'écrirai sur le contenu lié au traitement du signal acoustique en utilisant Python que je n'ai pas écrit dans article de l'année dernière.

table des matières

• Je souhaite convertir le taux d'échantillonnage
• Je souhaite lire un fichier audio 24 bits avec le module wave
• Je veux estimer le délai de groupe avec scipy.signal.group_delay
• Je veux encadrer en utilisant l'astuce de la foulée de numpy (je veux dessiner un spectrogramme)
• Je souhaite générer un signal TSP
• Je veux visualiser les zéros et les pôles

Je veux convertir le taux d'échantillonnage

Pour convertir la fréquence d'échantillonnage, si la fréquence d'échantillonnage d'origine et la fréquence d'échantillonnage que vous souhaitez convertir sont des ratios entiers, vous devez uniquement effectuer un échantillonnage haut / bas, mais dans le cas de ratios rationnels, utilisez le multiple commun minimum de chaque fréquence d'échantillonnage. Vous devez suréchantillonner puis sous-échantillonner à la fréquence d'échantillonnage souhaitée.

Nous envisageons ici de convertir le taux d'échantillonnage de 48 kHz à 44,1 kHz. La conversion du taux d'échantillonnage peut être réalisée avec le code suivant. La cible est les 10 premières secondes d'une source sonore à 48 kHz.

#!/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()

Puisque 44100/48000 = 147/160, la source sonore originale est d'abord suréchantillonnée jusqu'à un score de 147 fois. Ici, un tampon d'une taille de 147 fois est créé et le signal d'origine est remplacé dans le tampon à des intervalles de 147 points. Par conséquent, les points qui n'ont pas été attribués sont remplis de zéro. Ensuite, concevez le LPF pour n'inclure que les fréquences inférieures à la fréquence Nyquist 44100Hz et utilisez un filtre pour appliquer le LPF afin de limiter la bande. Le point ici est que la fréquence de Nyquist pendant le suréchantillonnage est de 48000/2 = 24000Hz, mais après le sous-échantillonnage, 44100/2 = 22050Hz devient la fréquence de Nyquist, de sorte que le filtre de limitation de bande appliqué pendant le suréchantillonnage est de 22050 ou moins. Il est possible d'utiliser un seul filtre qui devient. Après limitation de la bande, après correction du retard dû au LPF, un sous-échantillonnage est effectué à des intervalles de 160 points. En conséquence, une conversion de fréquence d'échantillonnage a été réalisée.

La forme d'onde audio et le spectrogramme réels sont les suivants. Les deux premiers sont après la conversion de 44,1 kHz et les deux derniers sont avant la conversion.

Je souhaite lire un fichier audio 24 bits avec le module wave

Si vous souhaitez lire un fichier audio 16 bits avec un module wave standard, vous pouvez le lire avec wave.open et en lire autant que vous le souhaitez avec des readframes et le convertir en int16 avec np.frombuffer, mais dans le cas d'une source sonore 24 bits, vous ne pouvez pas spécifier 24 bits avec frombuffer. Par conséquent, vous devez le lire vous-même. Ici, comme indiqué dans le code ci-dessous, lors de la lecture de 3 octets à la fois en utilisant unpack du module struct, la lecture d'une source sonore 24 bits est réalisée en compressant 0 et en décompressant en int32.

À propos, la source sonore en cours de lecture est source sonore d'exemple e-onkyo.

#!/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))

    #En emballant 0 dans le bit le plus bas et en déballant int
    #Extraire la valeur avec la valeur 24 bits en tant que entier 32 bits
    #  (<je suppose une petite valeur int andienne)
    #unpack renvoie un tuple[0]je prends le
    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.Normalisé à 0
    float_buf = np.where(ndarr_buf > 0,
                         ndarr_buf / (2.0 ** 31 - 1),
                         ndarr_buf / (2.0 ** 31))

    #Résoudre l'entrelacement(Pour source sonore stéréo)
    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()

Je souhaite estimer le délai de groupe avec scipy.signal.groupdelay

[Cet article] de l'article de l'année dernière (http://qiita.com/wrist/items/5759f894303e4364ebfd#%E7%BE%A4%E9%81%85%E5%BB%B6%E3%82%92%E8% Je n'ai pas remarqué quand j'ai écrit A8% 88% E7% AE% 97% E3% 81% 97% E3% 81% 9F% E3% 81% 84), mais je ne l'ai pas remarqué, mais la [méthode group_delay] pour trouver le retard de groupe ( https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.signal.group_delay.html#scipy.signal.group_delay) a été ajouté depuis la version 0.16.0 de scipy.signal. C'était. Cela peut être utilisé pour déterminer le retard de groupe d'un filtre numérique.

#!/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()

Je veux faire un traitement d'image en utilisant l'astuce de la foulée de numpy (je veux dessiner un spectrogramme)

En utilisant l'astuce de la foulée introduite dans les recettes 4.6 et 4.7 de IPython Data Science Cookbook, un tampon est divisé en images courtes. Vous pouvez effectuer efficacement un traitement d'image auquel on accède à plusieurs reprises sans générer une copie de l'image. En utilisant as_strided sous numpy.lib.stride_tricks, vous pouvez écrire la méthode de calcul du spectrogramme comme suit.

#!/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()

stride est une valeur qui représente l'intervalle d'octets lors de l'accès à une zone contiguë. Si vous modifiez cela, vous pourrez contrôler le nombre d'octets que vous avancez lorsque vous avancez d'un index de ligne ou de colonne.

as_strided est une méthode pour accéder au tampon avec cette foulée modifiée de manière pseudo. Le premier argument est le tampon de destination d'accès, le deuxième argument est la forme du tampon créé en modifiant la foulée et le troisième argument est la foulée que vous souhaitez modifier. Si vous accédez après avoir changé la mémoire tampon du premier argument à la foulée du troisième argument de manière pseudo, l'opération sera telle qu'un alias accessible en tant que tampon indiqué par la forme du deuxième argument est renvoyé.

Le format de cette foulée est (espacement des octets lors de l'avancement des lignes, espacement des octets lors de l'avancement des colonnes). Dans le cas du script ci-dessus, l'argument d'origine xs est un tableau unidimensionnel et sa foulée est (8,), donc si vous avancez l'index de 1, il avancera de 8 octets (dtype est float64). Si cette foulée est changée en (512 * 8, 8) par, par exemple, as_stride, xs, qui était à l'origine un tableau unidimensionnel, peut être considéré comme un tableau bidimensionnel xs_split. Si vous accédez à plusieurs reprises à xs_split dans une boucle for, chaque fois que vous y accédez, vous avancerez dans le sens de la ligne de xs_split, c'est-à-dire 512 * 8 octets avant xs, vous pouvez donc y accéder en décalant xs de 512pt. De cette manière, le traitement de trame qui serait autrement traité en découpant un tableau unidimensionnel en plusieurs trames de courte durée peut maintenant être effectué par des tours de pas sans générer une copie au moment de la découpe. Je vais.

Notez que cette astuce de foulée facilite l'accès en dehors de la plage du tampon interne, donc si vous faites une erreur, vous accéderez souvent à la zone non définie et tomberez.

Je souhaite générer un signal TSP

Lors de la conception d'un Sweep Pulse (signal TSP) utilisé pour la mesure acoustique, je pense qu'il est souvent obtenu en IFFTing ce qui a été conçu dans la région de fréquence et en le renvoyant dans la région de temps. Il peut être généré avec le code suivant. N_exp et m_exp, qui contribuent à la longueur effective, doivent être ajustés.

#!/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()

Le fichier wav généré par ceci est affiché comme audace comme indiqué ci-dessous.

Le signal TSP est décrit en détail dans le matériel de formation de base du professeur Kaneda pour la mesure de réponse impulsionnelle ci-dessous. http://www.asp.c.dendai.ac.jp/ASP/IRseminor2016.pdf

Je veux visualiser les zéros et les pôles

Vous pouvez convertir le coefficient de filtre temporel en zéro et en pôle avec tf2zpk, mais pour visualiser cela, vous pouvez dessiner un cercle en utilisant add_patch comme indiqué ci-dessous, puis tracer sans lignes et avec des marqueurs. .. J'ai remarqué après avoir collé la figure, mais comme il s'agit d'un filtre FIR, il n'y a pas de pôle, mais veuillez l'ignorer.

#!/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()

Autre

Je voulais aussi écrire le contenu suivant, mais je n'avais pas assez de temps, alors j'aimerais le rajouter ou l'écrire dans un autre article.

• Je veux estimer la réponse impulsionnelle des signaux d'entrée / sortie
• Conversion du filtre numérique au format SOS
• Je veux trouver l'enveloppe du signal en utilisant la transformation de Hilbert

Concernant la conversion Hilbert, je n'ai écrit que l'explication suivante, donc je la posterai pour le moment.

Afin de trouver l'enveloppe qui trace en douceur l'amplitude maximale du signal, il existe une méthode de multiplication du signal par la valeur absolue ou la valeur carrée par LPF, mais aussi la méthode pour la trouver à partir du signal d'analyse obtenu par conversion de Hilbert. A été

référence https://jp.mathworks.com/help/dsp/examples/envelope-detection.html