独立成分分析 Lab
在此 notebook 中，我们将使用独立成分分析方法从三个观察结果中提取信号，每个观察结果都包含不同的原始混音信号。这个问题与 ICA 视频中解释的问题一样。
数据集
首先看看手头的数据集。我们有三个 WAVE 文件，正如我们之前提到的，每个文件都是混音形式。如果你之前没有在 python 中处理过音频文件，没关系，它们实际上就是浮点数列表。
首先加载第一个音频文件 ICA_mix_1.wav [点击即可聆听该文件]：
import numpy as np
import wave
该wave模块为WAV声音格式提供了一个方便的界面。它不支持压缩/解压缩，但它支持单声道/立体声。
# read the wave file
mix_1_wave = wave.open('ICA_mix_1.wav', 'rb')
open返回的是一个Wave_read的实例，通过调用它的方法读取WAV文件的格式和数据
- getparams：一次性返回所有的WAV文件的格式信息，它返回的是一个组元(tuple)：声道数, 量化位数（byte单位）, 采样频率, 采样点数, 压缩类型, 压缩类型的描述。wave模块只支持非压缩的数据，因此可以忽略最后两个信息：
- getnchannels, getsampwidth, getframerate, getnframes等方法可以单独返回WAV文件的特定的信息。
- readframes：读取声音数据，传递一个参数指定需要读取的长度（以取样点为单位），readframes返回的是二进制数据（一大堆bytes)，在Python中用字符串表示二进制数据：
# 读取格式信息
# (nchannels, sampwidth, framerate, nframes, comptype, compname)
params = mix_1_wave.getparams()
nchannels, sampwidth, framerate, nframes = params[:4]
该文件只有一个声道（因此是单声道）。帧率是 44100，表示每秒声音由 44100 个整数组成（因为文件是常见的 PCM 16 位格式，所以是整数）。该文件总共有 264515 个整数/帧，因此时长为：264515/44100
我们从该 wave 文件中提取帧，这些帧将属于我们将运行 ICA 的数据集：
Wave_read.readframes(n)
最多可读取n帧音频，并以一串字节的形式进行读取和返回。
# Extract Raw Audio from Wav File
signal_1_raw = mix_1_wave.readframes(nframes)
# print(signal_1_raw)
signal_1 = np.frombuffer(signal_1_raw, 'Int16')
print(signal_1)
[ 879 1268 1460 ... 9204 9553 10006]
D:\Anaconda3\lib\site-packages\ipykernel_launcher.py:4: DeprecationWarning: Numeric-style type codes are deprecated and will result in an error in the future.
after removing the cwd from sys.path.
signal_1 现在是一个整数列表，表示第一个文件中包含的声音。
'length: ', len(signal_1) , 'first 100 elements: ',signal_1[:100]
('length: ',
264515,
'first 100 elements: ',
array([ 879, 1268, 1460, 1756, 1943, 2216, 2407, 2668, 2866,
3106, 3308, 3546, 3752, 3981, 4175, 4395, 4588, 4790,
4966, 5146, 5292, 5436, 5550, 5643, 5717, 5759, 5790,
5798, 5789, 5756, 5713, 5649, 5576, 5478, 5381, 5267,
5146, 4999, 4856, 4682, 4502, 4308, 4097, 3875, 3637,
3380, 3107, 2825, 2514, 2194, 1847, 1472, 1087, 671,
227, -219, -691, -1176, -1666, -2167, -2669, -3179, -3668,
-4170, -4643, -5116, -5559, -5985, -6380, -6765, -7105, -7422,
-7706, -7955, -8163, -8339, -8470, -8557, -8600, -8618, -8585,
-8524, -8425, -8298, -8129, -7947, -7720, -7475, -7205, -6916,
-6606, -6266, -5922, -5556, -5165, -4774, -4353, -3922, -3476,
-3021], dtype=int16))
如果将此数组绘制成线形图，我们将获得熟悉的波形：
import matplotlib.pyplot as plt
%matplotlib inline
fs = mix_1_wave.getframerate()
timing = np.linspace(0, len(signal_1)/fs, num=len(signal_1))
plt.figure(figsize=(12, 2))
plt.title("Recording 1")
plt.plot(timing, signal_1, c="#3ABFE7")
plt.ylim(-35000, 35000)
plt.show()
现在我们可以按照相同的方式加载另外两个 wave 文件 ICA_mix_2.wav 和 ICA_mix_3.wav
mix_2_wave = wave.open('ICA_mix_2.wav', 'rb')
params2 = mix_2_wave.getparams()
nchannels2, sampwidth2, framerate2, nframes2 = params2[:4]
signal_2_raw = mix_2_wave.readframes(nframes2)
signal_2 = np.frombuffer(signal_2_raw, 'Int16')
mix_3_wave = wave.open('ICA_mix_3.wav', 'rb')
params3 = mix_3_wave.getparams()
nchannels3, sampwidth3, framerate3, nframes3 = params3[:4]
signal_3_raw = mix_3_wave.readframes(nframes3)
signal_3 = np.frombuffer(signal_3_raw, 'Int16')
timing2 = np.arange(0, nframes2) * (1.0 / framerate2)
timing3 = np.linspace(0, nframes3 / framerate3, nframes3)
plt.figure(figsize=(12, 2))
plt.title("Recording 2")
plt.plot(timing2, signal_2)
plt.ylim(-35000, 35000)
plt.show()
plt.figure(figsize=(12, 2))
plt.title("Recording 3")
plt.plot(timing3, signal_3)
plt.ylim(-35000, 35000)
plt.show()
D:\Anaconda3\lib\site-packages\ipykernel_launcher.py:5: DeprecationWarning: Numeric-style type codes are deprecated and will result in an error in the future.
"""
D:\Anaconda3\lib\site-packages\ipykernel_launcher.py:11: DeprecationWarning: Numeric-style type codes are deprecated and will result in an error in the future.
# This is added back by InteractiveShellApp.init_path()
读取所有三个文件后，可以通过 zip 运算创建数据集。
- 通过将 signal_1、signal_2 和 signal_3 组合成一个列表创建数据集
X
X = list(zip(signal_1, signal_2, signal_3))
# Let's peak at what X looks like
X[:10]
[(879, 157, 879),
(1268, 899, 2584),
(1460, 1610, 4081),
(1756, 2325, 5682),
(1943, 3022, 7147),
(2216, 3718, 8689),
(2407, 4386, 10114),
(2668, 5056, 11583),
(2866, 5690, 12954),
(3106, 6318, 14327)]
现在准备运行 ICA 以尝试获取原始信号。
- 导入 sklearn 的 FastICA 模块
- 初始化 FastICA，查看三个成分
- 使用 fit_transform 对数据集 X 运行 FastICA 算法
# TODO: Import FastICA
from sklearn.decomposition import FastICA
# TODO: Initialize FastICA with n_components=3
fastICA = FastICA(n_components=3, random_state=0)
# TODO: Run the FastICA algorithm using fit_transform on dataset X
ica_result = fastICA.fit_transform(X)
ica_result.shape
(264515, 3)
new_result_signal_1 = ica_result[:, 0]
new_result_signal_2 = ica_result[:, 1]
new_result_signal_3 = ica_result[:, 2]
我们对信号进行绘制，查看波浪线的形状
plt.figure(figsize=(12, 15))
plt.subplot(311)
plt.title("Independent Component #1")
plt.plot(new_result_signal_1, c="#df8efd")
plt.ylim(-0.010, 0.010)
plt.subplot(312)
plt.title("Independent Component #2")
plt.plot(new_result_signal_2, c="#87de72")
plt.ylim(-0.010, 0.010)
plt.subplot(313)
plt.title("Independent Component #3")
plt.plot(new_result_signal_3, c="#f65e97")
plt.ylim(-0.010, 0.010)
plt.show()
某些波浪线看起来像音乐波形吗？
确认结果的最佳方式是聆听生成的文件。另存为 wave 文件并进行验证。在此之前，我们需要：
- 将它们转换为整数（以便另存为 PCM 16 位 Wave 文件），否则只有某些媒体播放器能够播放它们
- 将值映射到 int16 音频的相应范围内。该范围在 -32768 到 +32767 之间。基本的映射方法是乘以 32767。
- 音量有点低，我们可以乘以某个值（例如 100）来提高音量
from scipy.io import wavfile
# Convert to int, map the appropriate range, and increase the volume a little bit
result_signal_1_int = np.int16(new_result_signal_1*32767*100)
result_signal_2_int = np.int16(new_result_signal_2*32767*100)
result_signal_3_int = np.int16(new_result_signal_3*32767*100)
# Write wave files
wavfile.write("new_result_signal_1.wav", fs, result_signal_1_int)
wavfile.write("new_result_signal_2.wav", fs, result_signal_2_int)
wavfile.write("new_result_signal_3.wav", fs, result_signal_3_int)
现在生成的文件包括：[注意：确保先调低音响设备的音量，以防某些问题导致文件听起来像静态的]
音乐：
- 钢琴 - The Carnival of the Animals - XIII. The Swan（钢琴曲独奏版）。演奏者：Markus Staab
- 大提琴 - Cello Suite no. 3 in C, BWV 1009 - I. Prelude。演奏者：European Archive
扫一扫
1263
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
