女儿情：对于一段音乐进行频率频谱分析
【摘要】
简 介： 对于一段有筷子、饭盒、橡皮筋组成的乐器演奏的音乐进行分析。对于其音节的准确性讨论。可以看出，该乐器在高音所对应的频率节点还是比较准确，但低音中的（6,7,1）相对偏差比较大。令人惊奇的...
简 介： 对于一段有筷子、饭盒、橡皮筋组成的乐器演奏的音乐进行分析。对于其音节的准确性讨论。可以看出，该乐器在高音所对应的频率节点还是比较准确，但低音中的（6,7,1）相对偏差比较大。令人惊奇的是，居然伴音7,4频率点恰到好处，这究竟是什么魔法？
关键词： 频谱分析，音节
§01音乐频谱
1.1 音乐来源
今天在微信中友人给我发送了一段有趣的视频。这段视频演示了清新优雅的弹拨乐来自于饭盒上的九条橡皮筋，令人不禁感**:慨飞花摘叶皆可伤人，草木竹石均可为剑，诚不我欺也.**
▲ 图1.1.1 浪漫的音乐
▲ 图1.1.2 西游记主题曲：女儿情
剩下的问题就是对这段浪漫的音乐进行分析。
在 如何从MP4视频文件中抽取MP3音频？ 利用 格式工程 、 ffmpeg 软件从朋友圈下载的MP4视频中提取了相应的音频。由于原来的视频中音量比较小，使用 Audacity 软件对音频进行了初步处理，提升了音频数据音量。
1.2 读取和显示声音波形
1.2.1 读取波形代码
利用下面代码读取波形文件，并进行显示。
import sys,os,math,time
import matplotlib.pyplot as plt
from numpy import *
from scipy.io import wavfile
wavefile = '/home/aistudio/work/wave.wav'
sample_rate,sig = wavfile.read(wavefile)
print("sig.shape: {}".format(sig.shape), "sample_rate: {}".format(sample_rate))
plt.clf()
plt.figure(figsize=(10,6))
plt.plot(sig[:,0])
plt.xlabel("n")
plt.ylabel("wave")
plt.grid(True)
plt.tight_layout()
plt.savefig('/home/aistudio/stdout.jpg')
plt.show()
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
1.2.2 波形文件参数
可以看到波形文件的采样率以及双声道数据长度。
sig.shape: (1012608, 2)
sample_rate: 44100
- 1
- 2
1.3 波形时频联合分布
1.3.1 处理代码
from scipy import signal
f,t,Sxx = signal.spectrogram(sig[:, 0],
fs=sample_rate,
nperseg=8192,
noverlap=8000,
nfft=8192)
thread = 150000
Sxx[where(Sxx >thread)] = thread
startf = 25
endf = 200
plt.clf()
plt.figure(figsize=(15,7))
plt.pcolormesh(t, f[startf:endf], Sxx[startf:endf,:])
plt.xlabel('Time(s)')
plt.ylabel('Frequency(Hz)')
plt.grid(True)
plt.tight_layout()
plt.savefig('/home/aistudio/stdout.jpg')
plt.show()
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
1.3.2 显示处理结果
▲ 图1.3.1 音乐的是时频联合分布
thread = 50000
Sxx[where(Sxx >thread)] = thread
startf = 0
endf = 350
- 1
- 2
- 3
- 4
▲ 图1.3.2 音乐的是时频联合分布
▲ 图1.3.3 音频中的基频和各次谐波
▲ 声音前5秒钟的时频联合分布
▲ 图1.3.5 前两秒中的频谱
（1）完整的显示程序
from headm import * # =
from scipy.io import wavfile
wavefile = '/home/aistudio/work/wave.wav'
sample_rate,sig = wavfile.read(wavefile)
printt(sig.shape:, sample_rate:)
plt.clf()
plt.figure(figsize=(10,6))
plt.plot(sig[:,0])
plt.xlabel("n")
plt.ylabel("wave")
plt.grid(True)
plt.tight_layout()
plt.savefig('/home/aistudio/stdout.jpg')
plt.show()
from scipy import signal
f,t,Sxx = signal.spectrogram(sig[:sample_rate*5, 0],
fs=sample_rate,
nperseg=8192,
noverlap=8000,
nfft=8192)
thread = 50000
Sxx[where(Sxx >thread)] = thread
startf = 0
endf = 350
plt.clf()
plt.figure(figsize=(15,10))
plt.pcolormesh(t, f[startf:endf], Sxx[startf:endf,:])
plt.xlabel('Time(s)')
plt.ylabel('Frequency(Hz)')
plt.grid(True)
plt.tight_layout()
plt.savefig('/home/aistudio/stdout.jpg')
plt.show()
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
1.4 动态显示波形和频谱
1.4.1 音乐数据波形
下面是音乐的数据波形。可以看到背景的噪声的幅值还是比较大的。
▲ 图1.4.1 音乐数据波形
step_length = sample_rate//10
page_number = 300
overlap_length = step_length // 10
gifpath = '/home/aistudio/GIF'
if not os.path.isdir(gifpath):
os.makedirs(gifpath)
gifdim = os.listdir(gifpath)
for f in gifdim:
fn = os.path.join(gifpath, f)
if os.path.isfile(fn):
os.remove(fn)
from tqdm import tqdm
for id,i in tqdm(enumerate(range(page_number))):
startid = i * overlap_length
endid = startid + step_length
musicdata = sig[startid:endid,0]
plt.clf()
plt.figure(figsize=(10,6))
plt.plot(musicdata, label='Start:%d'%startid)
plt.xlabel("Samples")
plt.ylabel("Value")
plt.axis([0,step_length, -15000, 15000])
plt.grid(True)
plt.legend(loc='upper right')
plt.tight_layout()
savefile = os.path.join(gifpath, '%03d.jpg'%id)
plt.savefig(savefile)
plt.close()
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
1.4.2 音乐频谱变化
下面显示了音乐数据的频谱变化情况。
- 短时FFT参数：
-
显示频率：10~2000Hz
时间：100ms
补零长度：400ms
▲ 图1.4.2 音频的频谱变化数据
startf = 10
endf = 2000
startfid = int(step_length * startf / sample_rate*5)
endfid = int(step_length * endf / sample_rate*5)
from tqdm import tqdm
for id,i in tqdm(enumerate(range(page_number))):
startid = i * overlap_length
endid = startid + step_length
musicdata = list(sig[startid:endid,0])
zerodata = [0]*(step_length*4)
musicdata = musicdata + zerodata
mdfft = abs(fft.fft(musicdata))/step_length
fdim = linspace(startf, endf, endfid-startfid)
plt.clf()
plt.figure(figsize=(10,6))
plt.plot(fdim, mdfft[startfid:endfid], linewidth=3, label='Start:%d'%startid)
plt.xlabel("Frequency(Hz)")
plt.ylabel("Spectrum")
plt.axis([startf,endf, 0, 2000])
plt.grid(True)
plt.legend(loc='upper right')
plt.tight_layout()
savefile = os.path.join(gifpath, '%03d.jpg'%id)
plt.savefig(savefile)
plt.close()
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
§02 峰值频率
2.1 数据频率峰值频率
求出数据中频谱峰值频率的变化。
2.1.1 绘制所有频谱峰值
step_length = sample_rate//10
page_number = 500
overlap_length = step_length//5
startf = 10
endf = 2000
startfid = int(step_length * startf / sample_rate*5)
endfid = int(step_length * endf / sample_rate*5)
maxfdim = []
maxadim = []
maxtdim = []
for id,i in tqdm(enumerate(range(page_number))):
startid = i * overlap_length
endid = startid + step_length
if endid > sig.shape[0]: break
musicdata = list(sig[startid:endid,0])
zerodata = [0]*(step_length*4)
musicdata = musicdata + zerodata
mdfft = abs(fft.fft(musicdata))/step_length
fdim = linspace(startf, endf, endfid-startfid)
videofft = mdfft[startfid:endfid]
maxid = argmax(videofft)
maxfreq = fdim[maxid]
maxampl = videofft[maxid]
starttime = startid / sample_rate
maxfdim.append(maxfreq)
maxadim.append(maxampl)
maxtdim.append(starttime)
col1 = 'steelblue'
col2 = 'red'
plt.figure(figsize=(12,6))
plt.plot(maxtdim, maxfdim, color=col1, linewidth=3)
plt.xlabel('Time', fontsize=14)
plt.ylabel('Max Frequency(Hz)', color=col1, fontsize=14)
plt.grid(True)
ax2 = plt.twinx()
ax2.plot(maxtdim, maxadim, color=col2, linewidth=1)
ax2.set_ylabel('Pick Value', color=col2, fontsize=16)
plt.tight_layout()
plt.savefig('/home/aistudio/stdout.jpg')
plt.show()
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
下图给出了频谱中峰值频率以及峰值幅度的变化
▲ 图2.1.1 音乐数据频率峰值
▲ 图2.1.2 音乐数据频率峰值
▲ 图2.1.3 音乐数据频率峰值
2.1.2 浏览频谱
▲ 图2.1.4 数据频谱峰值频率与幅值变化
from headm import * # =
datafile = '/home/aistudio/maxdata.npz'
data = load(datafile)
printt(data.files)
maxfdim = data['maxf']
maxadim = data['maxa']
maxtdim = data['maxtdim']
gifpath = '/home/aistudio/GIF'
if not os.path.isdir(gifpath):
os.makedirs(gifpath)
gifdim = os.listdir(gifpath)
for f in gifdim:
fn = os.path.join(gifpath, f)
if os.path.isfile(fn):
os.remove(fn)
page_number = 200
data_length = len(maxfdim)
step_length = data_length // 10
for i in tqdm(range(page_number)):
startid = int(i * (data_length - step_length) / page_number)
endid = startid + step_length
if endid >= data_length: endid = data_length
plt.clf()
col1 = 'steelblue'
col2 = 'red'
plt.figure(figsize=(12,8))
plt.plot(maxtdim[startid:endid], maxfdim[startid:endid], color=col1, linewidth=3)
plt.axis([maxtdim[startid], maxtdim[endid-1], 0, 2000])
plt.xlabel('Time', fontsize=14)
plt.ylabel('Max Frequency(Hz)', color=col1, fontsize=14)
plt.grid(True)
ax2 = plt.twinx()
ax2.plot(maxtdim[startid:endid], maxadim[startid:endid], color=col2, linewidth=1)
plt.axis([maxtdim[startid], maxtdim[endid-1], 0, 5000])
ax2.set_ylabel('Pick Value', color=col2, fontsize=16)
plt.tight_layout()
savefile = os.path.join(gifpath, '%03d.jpg'%i)
plt.savefig(savefile)
plt.close()
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
2.2 频率直方图
下面给出了在音乐中的频率峰值的统计。
plt.clf()
plt.figure(figsize=(12,8))
plt.hist(maxfdim, bins=500)
plt.xlabel('Frequency')
plt.ylabel('Histogram')
plt.grid(True)
plt.savefig('/home/aistudio/stdout.jpg')
plt.show()
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
▲ 图2.2.1 频率峰值的统计
▲ 在 350Hz ~ 1000Hz 之间的直方图统计
▲ 图2.2.3 根据乐谱获得频率点
2.3 所有音节频率
以 音符3 所对应的频率 f 3 = 695.9 H z f_3 = 695.9Hz f3=695.9Hz ，为标准，计算所有音节的位置频率：
bw = where((bins > 350) & (bins < 1000))
fdim = frequency[bw]
bdim = bins[bw]
fmaxid = argmax(fdim)
print(bdim[fmaxid])
- 1
- 2
- 3
- 4
- 5
step12 = [-9, -7, -5, -4,-2,0,1,3,5]
stepnum = [5,6,7,1,2,3,4,5,6]
freq3 = 695.9
for id,s in enumerate(step12):
freq = 2**(s/12)*freq3
print('%d : %.2fHz'%(stepnum[id], freq))
- 1
- 2
- 3
- 4
- 5
- 6
- 7
5 : 413.78Hz
6 : 464.46Hz
7 : 521.34Hz
1 : 552.34Hz
2 : 619.98Hz
3 : 695.90Hz
4 : 737.28Hz
5 : 827.57Hz
6 : 928.92Hz
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
▲ 图2.3.1 音节频率值
step12 = [-9, -7, -5, -4,-2,0,1,3,5]
stepnum = [5,6,7,1,2,3,4,5,6]
freq3 = 695.9
plt.clf()
plt.figure(figsize=(12,8))
for id,s in enumerate(step12):
freq = 2**(s/12)*freq3
print('%d : %.2fHz'%(stepnum[id], freq))
plt.plot([freq,freq],[0, 250], linewidth=2, c='red')
plt.plot(bins[bw], frequency[bw])
plt.xlabel('Frequency')
plt.ylabel('Histogram')
plt.grid(True)
plt.savefig('/home/aistudio/stdout.jpg')
plt.show()
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
▲ 图2.3.2 根据乐谱获得频率点
※ 总 结 ※
对于一段有筷子、饭盒、橡皮筋组成的乐器演奏的音乐进行分析。对于其音节的准确性讨论。可以看出，该乐器在高音所对应的频率节点还是比较准确，但低音中的（6,7,1）相对偏差比较大。
令人惊奇的是，居然伴音7,4频率点恰到好处，这究竟是什么魔法？
■ 相关文献链接:
● 相关图表链接:
- 图1.1.1 浪漫的音乐
- 图1.1.2 西游记主题曲：女儿情
- 图1.3.1 音乐的是时频联合分布
- 图1.3.2 音乐的是时频联合分布
- 图1.3.3 音频中的基频和各次谐波
- 声音前5秒钟的时频联合分布
- 图1.3.5 前两秒中的频谱
- 图1.4.1 音乐数据波形
- 图1.4.2 音频的频谱变化数据
- 图2.1.1 音乐数据频率峰值
- 图2.1.2 音乐数据频率峰值
- 图2.1.3 音乐数据频率峰值
- 图2.1.4 数据频谱峰值频率与幅值变化
- 图2.2.1 频率峰值的统计
- 在 350Hz ~ 1000Hz 之间的直方图统计
- 图2.2.3 根据乐谱获得频率点
- 图2.3.1 音节频率值
- 图2.3.2 根据乐谱获得频率点
文章来源: zhuoqing.blog.csdn.net，作者：卓晴，版权归原作者所有，如需转载，请联系作者。
原文链接：zhuoqing.blog.csdn.net/article/details/122912698
【版权声明】本文为华为云社区用户转载文章，如果您发现本社区中有涉嫌抄袭的内容，欢迎发送邮件进行举报，并提供相关证据，一经查实，本社区将立刻删除涉嫌侵权内容，举报邮箱：
cloudbbs@huaweicloud.com
- 点赞
- 收藏
- 关注作者
评论（0）