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

pythonでフーリエ変換をしてみる

pythonでいろいろ計算処理をできるようになるため、まずはフーリエ変換してグラフに表示するところから始めてみる。

検索すれば山ほどフーリエ変案の例がでてくるが実際に自分で手を動かしてやってみる。

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt


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

    return tmp

    
# 信号にノイズを加える
def add_noise(sig, sn):
    # 信号の電力(dB)
    norm_sig=np.linalg.norm(sig)
    pw_sig_dB=20*np.math.log10(norm_sig)
    # ノイズ生成
    sig_length=len(sig)
    noise=np.random.randn(sig_length)
    pw_noise_dB=20*np.math.log10(np.linalg.norm(noise))
    # ノイズの電力を指定のSNになるように調整する
    pw_noise_dB_delta=(pw_sig_dB-sn)-pw_noise_dB
    noise=noise*10**(pw_noise_dB_delta/20)
    sig_with_noise=sig+noise

    return sig_with_noise

if __name__ == '__main__':
    sa=16000
    fftpt=2048
    sn=10
    sig=np.exp(2*np.pi*1j*(100/fftpt)*np.arange(0,3000))
    norm_sig=np.linalg.norm(sig)
    sig=sig/norm_sig
    sig = add_noise(sig, sn)
    
    #フーリエ変換
    w=np.kaiser(fftpt,0.5)
    fftdt=np.fft.fft(w*sig[:fftpt])/np.sqrt(fftpt) #matlabでいうfftshiftがかかった状態
    fftdt_pw=20*np.log10(np.abs(fftdt))
    freq_list = np.arange(0,fftpt)*sa/fftpt
    
    # 表示のためにノイズレベル、最大値を求める
    # 最頻値をノイズレベルと仮定する
    sudo_noise_level_dB = np.median(fftdt_pw) 
    sudo_noise_level_dB =ciel_quantize(sudo_noise_level_dB, 10)
    # 電力の最大値
    max_pw_dB = np.max(fftdt_pw)
    max_pw_dB = ciel_quantize(max_pw_dB, 10)

    #表示
    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]')
    plt.axis([freq_list[0], freq_list[-1], sudo_noise_level_dB, max_pw_dB])
    plt.show()

関連記事

コメント

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

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

スポンサード リンク

カテゴリー

スポンサード リンク