質問編集履歴

1

コード情報の追加

2020/06/09 04:13

投稿

退会済みユーザー
test CHANGED
File without changes
test CHANGED
@@ -14,8 +14,176 @@
14
14
 
15
15
  データは22050あります。
16
16
 
17
+ 現時点では、Excel1行目にA~XFD列まで数値が並ぶ状態です。
18
+
19
+
20
+
17
21
  ### 当該のソースコード
18
22
 
23
+ ```ここに言語を入力
24
+
25
+ import wave
26
+
27
+ import struct
28
+
29
+ from scipy import fromstring, int32
30
+
31
+ import numpy as np
32
+
33
+ from pylab import *
34
+
35
+ %matplotlib inline
36
+
37
+
38
+
39
+ wavfile = '440Hz.wav' #サンプルwavファイル
40
+
41
+ wr = wave.open(wavfile, "rb")
42
+
43
+ ch = wr.getnchannels() # モノラルなら1,ステレオなら2
44
+
45
+ width = wr.getsampwidth() # バイト数 (1byte=8bit)
46
+
47
+ fr = wr.getframerate() #サンプリンググレート(サンプリング周波数)
48
+
49
+ fn = wr.getnframes() # 総フレーム数(データ分割数)⇒サンプリング周波数で割れば時間
50
+
51
+
52
+
53
+ #wavfileから(N*span)のデータを先頭から切り出す
54
+
55
+ #連続的なデータを離散的に表すためには、データの成分は"fr"の1/2未満(ナイキスト周波数)でなければならない。
56
+
57
+ N = 22050 #サンプリングレート"fr"の半分の値
58
+
59
+ span =4 #Nの整数倍 #周波数分解能:分解できる周波数の細かさ
60
+
61
+
62
+
63
+ print('サンプル数',N)
64
+
65
+ print('チャンネル', ch)
66
+
67
+ print('バイト数', width)
68
+
69
+ print('サンプリンググレート', fr)
70
+
71
+ print('総フレーム数', fn)
72
+
73
+ print('時間',fn/fr,'秒')
74
+
75
+ #"N*span"分の時間の取得
76
+
77
+ print('サンプル時間', 1.0 * N * span / fr, '秒')
78
+
79
+
80
+
81
+ origin = wr.readframes(wr.getnframes()) #メソッドreadframes(n)でnフレーム分のデータを読み込む、ここでは総フレーム数分の読み込み
82
+
83
+ data = origin[:N * span * ch * width] #"origin"から要素の範囲[ ]を指定
84
+
85
+ wr.close()
86
+
87
+
88
+
89
+ print('現配列長', len(origin)) #"origin"の要素数
90
+
91
+ print('サンプル配列長: ', len(data)) #"data"の要素数
92
+
93
+ ```
94
+
95
+ サンプル数 22050
96
+
97
+ チャンネル 1
98
+
99
+ バイト数 2
100
+
101
+ サンプリンググレート 44100
102
+
103
+ 総フレーム数 441000
104
+
105
+ 時間 10.0 秒
106
+
107
+ サンプル時間 2.0 秒
108
+
109
+ 現配列長 882000
110
+
111
+ サンプル配列長: 176400
112
+
113
+ ```ここに言語を入力
114
+
115
+ X = np.frombuffer(data, dtype="int16")#"data"をバイナリ表記から16bitsの整数数列に変換
116
+
117
+
118
+
119
+ # ステレオ前提、左右音に分ける ※モノラルは単に1つおきにデータを読みこむため、必要ない工程
120
+
121
+ left = X[::2] #"0から2番目おき"に要素を得る
122
+
123
+ right = X[1::2] #"1から2番目おき"に要素を得る
124
+
125
+ ```
126
+
127
+ ```ここに言語を入力
128
+
129
+ def fourier (x, n, w): #x:データ成分、n:個数、w:次元
130
+
131
+ K = []
132
+
133
+ for i in range(0, w-2): #2番目おきに要素を得ているため、0~N-2範囲となる
134
+
135
+ sample = x[i * n:( i + 1) * n] #i~(i+1)番目の要素を得る
136
+
137
+ partial = np.fft.fft(sample) #"sample"をフーリエ変換
138
+
139
+ K.append(partial) #"K"に"partial"を追加
140
+
141
+
142
+
143
+ return K
144
+
145
+
146
+
147
+ #周波数分布をもとに、実空間での波形を生成しています
148
+
149
+ def inverse_fourier (k):
150
+
151
+ ret = []
152
+
153
+ for sample in k:
154
+
155
+ inv = np.fft.ifft(sample) #"sample"を逆フーリエ変換
156
+
157
+ ret.extend(inv.real) #"inv.real"を"ret"に追加
158
+
159
+
160
+
161
+ print (len(sample))
162
+
163
+ return ret
164
+
165
+
166
+
167
+ Kl = fourier(left, N, span)
168
+
169
+ Kr = fourier(right, N, span)
170
+
171
+ freqlist = np.fft.fftfreq(N, d=1/fr) #周波数リスト
172
+
173
+ amp = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in Kl[1]] #振幅スペクトル #実部と虚部を取り出すには、".real" と ".imag" を使用
174
+
175
+ plot(freqlist, amp, marker= 'o', linestyle='-') #周波数リスト、振幅スペクトル、点、線スタイル
176
+
177
+ axis([0, 1000 , 0, 100000])
178
+
179
+
180
+
181
+ amp = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in Kr[1]]
182
+
183
+ plot(freqlist, amp, marker= 'o', linestyle='-')
184
+
185
+ ```
186
+
19
187
  ```
20
188
 
21
189
  amp = [[np.sqrt(c.real ** 2 + c.imag ** 2) for c in Kl[1]]]