いろんなものはつながっている

pythonでフィルタをかけてみる

16kサンプリングのwavデータを読み込んで、8kに帯域をカットして、8kサンプリングのwavデータを保存する。

# -*- coding: utf-8 -*-
import os.path
import struct
import wave
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt


# WAVファイルを読み込む
def read_wav(file_path):
    with wave.open(file_path , "rb" ) as wf:
        #チャンネル数(モノラル/ステレオ)
        channel_num = wf.getnchannels()
        #サンプルサイズ(バイト)
        sample_size = wf.getsampwidth()
        #サンプリング周波数(Hz)
        sampling_rate = wf.getframerate()
        #フレーム数
        frame_num = wf.getnframes()

        #データ読み込み
        wav_date = wf.readframes(frame_num)
        wav_date_int16 = np.frombuffer(wav_date, dtype= "int16")
        wav_date_float = wav_date_int16 / 32767.0

        return wav_date_float, sampling_rate


# WAVファイル書き込み
def write_wav(file_path, wav_data, sampling_rate, 
              sample_size = 2, channel_num = 1):
    with wave.open(file_path, "wb" ) as wr:
        wr.setnchannels(channel_num)
        wr.setsampwidth(sample_size)
        wr.setframerate(sampling_rate)
        wr.writeframes(wav_data)


# 指定桁で切り上げする
def ciel_quantize(src, digit):
    tmp = src/digit;
    tmp = np.math.trunc(tmp)
    if tmp >= 0: 
        tmp = (int)((tmp + 1) * digit)
    else:
        tmp = (int)((tmp - 1) * digit)

    return tmp    


# フーリエ変換
def fft_power(sig, fftpt):
    w=np.kaiser(fftpt,0.5)
    fftdt=np.fft.fft(w*sig[:fftpt])
    fftdt_pw=20*np.log10(np.abs(fftdt))

    return fftdt_pw

# フーリエ変換結果表示
def show_fft(fftdt_pw, sa, fftpt, show = True, complex_data = True, yaxis_min_type = 'noise_level'):
    # Y軸表示の最小値、最大値を求める
    min_pw_dB = 0
    if yaxis_min_type == 'noise_level':
        sudo_noise_level_dB = np.median(fftdt_pw) 
        sudo_noise_level_dB =ciel_quantize(sudo_noise_level_dB, 10)
        min_pw_dB = sudo_noise_level_dB
    else :
        min_pw_dB = np.min(fftdt_pw)
        min_pw_dB = ciel_quantize(min_pw_dB, 10)
            
    # 最頻値をノイズレベルと仮定する
    # 電力の最大値
    max_pw_dB = np.max(fftdt_pw)
    max_pw_dB = ciel_quantize(max_pw_dB, 10)

    #周波数軸
    freq_list = np.arange(0,fftpt) * (sa/fftpt)

    #表示
    # plt.title('FFT')
    plt.plot(freq_list, fftdt_pw)
    fft_freq_resolution = '%.2f' % (sa/fftpt)
    plt.xlabel('frequency [Hz]  ' + fft_freq_resolution + 'Hz/bin')
    plt.ylabel('power[dB]')
    freq_list_start = freq_list[0]
    freq_list_end = 0
    if complex_data:
        freq_list_end = freq_list[-1]
    else:
        freq_list_end = freq_list[int(len(freq_list)/2)]
    plt.axis([freq_list_start, freq_list_end, min_pw_dB, max_pw_dB])

    if show:
        plt.show()


if __name__ == '__main__':
    # WAVファイル読み込み
    wave_file_path=r'C:\audiodata_16k.wav'
    wave_date, sampling_rate = read_wav(wave_file_path)

    # フーリエ変換
    fftpt=2048
    fftdt = fft_power(wave_date, fftpt)

    # フーリエ変換結果表示
    plt.subplot(3, 1, 1)
    show_fft(fftdt, sampling_rate, fftpt, False, False, 'min')

    # フィルタ
    # 参考:https://org-technology.com/posts/low-pass-filter.html
    numtaps=255
    cutoff = 1.0 / 2.0 #帯域を16k->8kへ
    b = signal.firwin(numtaps, cutoff) 
    wave_date_filtered = signal.lfilter(b, 1, wave_date)
    wave_date_filtered = wave_date_filtered[numtaps:]

    # フーリエ変換
    fftdt_filtered = fft_power(wave_date_filtered, fftpt)

    # フーリエ変換結果表示
    plt.subplot(3, 1, 2)
    show_fft(fftdt_filtered, sampling_rate, fftpt, False, False, 'min')

    # サンプリング周波数を半分にする
    wav_data_down_sample = wave_date_filtered[::2]

    # フーリエ変換
    fftdt_down_sampling = fft_power(wav_data_down_sample, fftpt)

    plt.subplot(3, 1, 3)
    show_fft(fftdt_down_sampling, sampling_rate/2, fftpt, False, False, 'min')
    plt.tight_layout()

    # wavファイルに出力する
    # そのためにバイナリに変換する
    wave_data_ = [np.int16(x * 32767) for x in wav_data_down_sample]
    wave_data_ = struct.pack("h" * len(wave_data_), *wave_data_)
    root, ext = os.path.splitext(wave_file_path)
    file_path_created_wav = root + '_8k.wav'
    write_wav(file_path_created_wav, wave_data_, int(sampling_rate/2))
    

scipyを用いるとwavデータの読み込みをmatlabのようにできる。

#wavファイルの読み込み
from scipy.io.wavfile import read as audioread
fs, data = audioread(wave_file_path)

#wavファイルの書き込み
from scipy.io.wavfile import write as audiowrite
audiowrite(wave_file_path, fs, data)

関連記事

コメント

  1. この記事へのコメントはありません。

  1. この記事へのトラックバックはありません。

スポンサード リンク

カテゴリー

スポンサード リンク