2020-01-18 python實(shí)現(xiàn)stft并繪制時(shí)頻譜

官方文檔中給出了非常詳細(xì)的安裝方法

公司主營(yíng)業(yè)務(wù)：成都網(wǎng)站設(shè)計(jì)、網(wǎng)站制作、外貿(mào)營(yíng)銷網(wǎng)站建設(shè)、移動(dòng)網(wǎng)站開發(fā)等業(yè)務(wù)。幫助企業(yè)客戶真正實(shí)現(xiàn)互聯(lián)網(wǎng)宣傳，提高企業(yè)的競(jìng)爭(zhēng)能力。成都創(chuàng)新互聯(lián)公司是一支青春激揚(yáng)、勤奮敬業(yè)、活力青春激揚(yáng)、勤奮敬業(yè)、活力澎湃、和諧高效的團(tuán)隊(duì)。公司秉承以“開放、自由、嚴(yán)謹(jǐn)、自律”為核心的企業(yè)文化，感謝他們對(duì)我們的高要求，感謝他們從不同領(lǐng)域給我們帶來的挑戰(zhàn)，讓我們激情的團(tuán)隊(duì)有機(jī)會(huì)用頭腦與智慧不斷的給客戶帶來驚喜。成都創(chuàng)新互聯(lián)公司推出白河免費(fèi)做網(wǎng)站回饋大家。

函數(shù)聲明：

librosa.core.stft(y, n_fft=2048, hop_length=None, win_length=None, window='hann', center=True, dtype=class 'numpy.complex64', pad_mode='reflect')

常用參數(shù)說明：

y：輸入的numpy數(shù)組，要求都是實(shí)數(shù)

n_fft：fft的長(zhǎng)度，默認(rèn)2048

hop_length：stft中窗函數(shù)每次步進(jìn)的單位

win_length：窗函數(shù)的長(zhǎng)度

window：窗函數(shù)的類型

return：一個(gè)1+n_fft/2*1+len(y)/hop_length的二維復(fù)數(shù)矩陣，其實(shí)就是時(shí)頻譜

參考：

主要用這兩個(gè)

matplotlib.pyplot.pcolormesh()

matplotlib.pyplot.colorbar()

Python科學(xué)計(jì)算——復(fù)雜信號(hào)FFT

FFT (Fast Fourier Transform, 快速傅里葉變換) 是離散傅里葉變換的快速算法，也是數(shù)字信號(hào)處理技術(shù)中經(jīng)常會(huì)提到的一個(gè)概念。用快速傅里葉變換能將時(shí)域的數(shù)字信號(hào)轉(zhuǎn)換為頻域信號(hào)，轉(zhuǎn)換為頻域信號(hào)后我們可以很方便地分析出信號(hào)的頻率成分。

當(dāng)我們把雙頻信號(hào)FFT示例中的 fft_size 的值改為 2**12 時(shí)，這時(shí)，基頻為 16Hz，不能被 1kHz整除，所以 1kHz 處發(fā)生了頻譜泄露，而它能被 4kHz 整除，所以 4kHz 可以很好地被采樣。

由于波形的前后不是連續(xù)的，出現(xiàn)波形跳變，而跳變處有著非常廣泛的頻譜，因此FFT的結(jié)果中出現(xiàn)了頻譜泄漏。

為了減小FFT所截取的數(shù)據(jù)段前后的跳變，可以對(duì)數(shù)據(jù)先乘以一個(gè)窗函數(shù)，使得其前后數(shù)據(jù)能平滑過渡。常用的hanning窗函數(shù)的定義如下：

50Hz 正弦波與hann窗函數(shù)乘積之后的重復(fù)波形如下：

我們對(duì)頻譜泄漏示例中的1kHz 和 4kHz 信號(hào)進(jìn)行了 hann 窗函數(shù)處理，可以看出能量更加集中在 1kHz 和 4kHz，在一定程度上抑制了頻譜泄漏。

以 1kHz 三角波為例，我們知道三角波信號(hào)中含有豐富的頻率信息，它的傅里葉級(jí)數(shù)展開為：

當(dāng)數(shù)字信號(hào)的頻率隨時(shí)間變化時(shí)，我們稱之為掃頻信號(hào)。以頻率隨時(shí)間線性變化的掃頻信號(hào)為例，其數(shù)學(xué)形式如下：

其頻率隨時(shí)間線性變化，當(dāng)我們?cè)?[0,1] 的時(shí)間窗口對(duì)其進(jìn)行采樣時(shí)，其頻率范圍為 0~5kHz。當(dāng)時(shí)間是連續(xù)時(shí)，掃頻信號(hào)的頻率也是連續(xù)的。但是在實(shí)際的處理中，是離散的點(diǎn)采樣，因此時(shí)間是不連續(xù)的，這就使掃頻信號(hào)的快速傅里葉變換問題退化為多點(diǎn)頻信號(hào)快速傅里葉變換問題。其快速傅里葉變換得到的頻譜圖如下所示：

以 50Hz 正弦信號(hào)相位調(diào)制到 1kHz 的信號(hào)為例，其信號(hào)形式如下：

它的時(shí)域波形，頻率響應(yīng)和相位響應(yīng)如下圖所示：

以掃頻信號(hào)為例，當(dāng)我們要探究FFT中的能量守恒時(shí)，我們要回歸到信號(hào)最初的形式：

python 二維FFT

二維FFT常用在圖像處理上，首先要能理解二維FFT的意義，否則很難明白它到底是怎么工作的。

第一列是原圖和對(duì)應(yīng)的頻率信息，第二列是去除低頻部分后，F(xiàn)FT逆變換得到的圖像。第三列是去除高頻部分后FFT逆變換得到的圖像。

從第二列可以看出高頻貢獻(xiàn)了圖像的細(xì)節(jié)。從白到黑的邊界保留了下來。而原圖中大片的白與大片的黑在這個(gè)圖中沒什么區(qū)別。

第三列中保留了原圖中的亮部與灰部，而由黑到白的臨界線卻很模糊。細(xì)小的白線黑線也沒能顯示。所以低頻貢獻(xiàn)了圖像的明暗。

2.工作原理理解

二維FFT就是先對(duì)行做次一維FFT，這樣每個(gè)元素都是關(guān)于行頻率信息了，然后再對(duì)列做一維FFT,這樣每個(gè)元素都包含了行和列的頻率信息。每個(gè)元素都是個(gè)復(fù)數(shù)，取絕對(duì)值可得到振幅，從實(shí)部與虛部的比值可等到相位，在二維矩陣的位置信息包含了頻率大小和方向。方向在一維FFT中是不用考慮的。

FFT2的結(jié)果也是正頻率從0到高然后負(fù)頻率從高到0.fftshift()之后會(huì)將低頻放到中間位置。

第一幅圖的頻譜是中間一條白線，也就是說許多個(gè)正弦波沿橫向傳播。縱向上沒有變化。

第三幅圖的頻譜是十字形加一條從左下角到右上角的直線。說明原圖在橫向，縱向都有變化，變化的方向從左下角到右上角。

從中心到頻譜圖上某一點(diǎn)構(gòu)成的向量方向就是這個(gè)波傳播的方向。

正負(fù)對(duì)稱才能消除虛部，這點(diǎn)與一維FFT原理一致。

如何用python實(shí)現(xiàn)音樂的頻譜顯示

用pymedia可以很容易的實(shí)現(xiàn)：

import pymedia.audio.acodec as acodec

sName='YOUR FILENAME'

dec= acodec.Decoder( str.split( sName, '.' )[ -1 ].lower() )

f= open( sName, 'rb' )

s= f.read( 8192 )

r= dec.decode( s )

import pymedia.audio.sound as sound

snd= sound.Output( r.sample_rate, r.channels, sound.AFMT_S16_LE )

while len( s )0:

if r: snd.play( r.data )

s= f.read( 512 )

r= dec.decode( s )

import time

while snd.isPlaying(): time.sleep( .05 )

誰知道如何在python中用處理wav文件，并且對(duì)他的頻譜進(jìn)行分析的程序

1.讀取wav文件

# -*- coding: utf-8 -*-

import wave

import pylab as pl

import numpy as np

# 打開WAV文檔

f = wave.open(r"c:\WINDOWS\Media\ding.wav", "rb")

# 讀取格式信息

# (nchannels, sampwidth, framerate, nframes, comptype, compname)

params = f.getparams()

nchannels, sampwidth, framerate, nframes = params[:4]

# 讀取波形數(shù)據(jù)

str_data = f.readframes(nframes)

f.close()

#將波形數(shù)據(jù)轉(zhuǎn)換為數(shù)組

wave_data = np.fromstring(str_data, dtype=np.short)

wave_data.shape = -1, 2

wave_data = wave_data.T

time = np.arange(0, nframes) * (1.0 / framerate)

# 繪制波形

pl.subplot(211)

pl.plot(time, wave_data[0])

pl.subplot(212)

pl.plot(time, wave_data[1], c="g")

pl.xlabel("time (seconds)")

pl.show()

2.觀察信號(hào)頻譜

# -*- coding: utf-8 -*-

import numpy as np

import pylab as pl

sampling_rate = 8000

fft_size = 512

t = np.arange(0, 1.0, 1.0/sampling_rate)

x = np.sin(2*np.pi*156.25*t) + 2*np.sin(2*np.pi*234.375*t)

xs = x[:fft_size]

xf = np.fft.rfft(xs)/fft_size

freqs = np.linspace(0, sampling_rate/2, fft_size/2+1)

xfp = 20*np.log10(np.clip(np.abs(xf), 1e-20, 1e100))

pl.figure(figsize=(8,4))

pl.subplot(211)

pl.plot(t[:fft_size], xs)

pl.xlabel(u"時(shí)間(秒)")

pl.title(u"156.25Hz和234.375Hz的波形和頻譜")

pl.subplot(212)

pl.plot(freqs, xfp)

pl.xlabel(u"頻率(Hz)")

pl.subplots_adjust(hspace=0.4)

pl.show()