第六章 離散餘弦轉換 | Discrete Cosine Transform
這一章的主題是離散徐弦轉換 Discrete Cosine Transform (DCT)，這用在 MP3 與類似的壓縮音樂，JPEG 與類似的壓縮圖片，MPEG 類的影像格式。DCT 在許多方面都與 DFT (Discrete Fourier Transform) 相似，之前我們有用來做一些分析。一旦我們學會 DCT 的原理，也就比較容易解釋 DFT。
我們依循以下步驟來解釋：
- 我們從合成問題開始：給定一組頻率組成，與它們的振幅，我們該如何組成一個波？
- 接下來，我們重寫合成問題的程式，改用 NumPy array。這一步會增加效能，也同時為下一步提供透視。
- 我們會看分析問題：給定一個訊號與一組頻率，我們要如何找到每個頻率的振幅？我們會從一個概念比較簡單但比較慢的解法開始。
- 最後，我們會用幾個線性代數的概念來找到更有效率的演算法。如果你已經熟悉線性代數，那很棒。總之我會解釋所需要的事情。
6.1 合成 | Synthesis
假設我給你一串振幅與一串頻率，然後要你組成一個訊號，它是這些頻率成份的總和。可使用 thinkdsp 模組裡的物件，這有一個簡單的方法執行這個操作，就叫做合成 synthesis：def synthesize1(amps, fs, ts):
components = [thinkdsp.CosSignal(freq, amp)
for amp, freq in zip(amps, fs)]
signal = thinkdsp.SumSignal(*components)
ys = signal.evaluate(ts)
return yscomponents 是一串 CosSignal 物件，每一個振幅與頻率的對應配對。SumSignal 代表這些頻率成份的總和。
最後，evaluate 計算每一個在 ts 裡的時間所對應的訊號值。
我們可以測一下這個函式：
amps = np.array([0.6, 0.25, 0.1, 0.05])
fs = [100, 200, 300, 400]
framerate = 11025
ts = np.linspace(0, 1, framerate)
ys = synthesize1(amps, fs, ts)
wave = thinkdsp.Wave(ys, framerate)6.2 使用陣列合成 | Synthesis with arrays
圖6.1 用陣列合成
這裡有另個方式寫出 synthesize：
def synthesize2(amps, fs, ts):
args = np.outer(ts, fs)
M = np.cos(PI2 * args)
ys = np.dot(M, amps)
return ys- np.outer 計算 ts 與 fs 的外積(outer product)，結果是一個陣列，有一個橫列的對應 ts 的元素，有一個直行對應 fs 的元素。每一個元素在陣列裡就是頻率與時間的乘積，。
- 我們將 args 乘上 然後做 運算，所以結果的每個元素是 。因為 ts 是沿著直行下來，每個直行都包含了在某特定頻率的 cosine 訊號，其由時間演算而來。
- np.dot 把 M 的每個橫列與 amps 相乘(以元素)，然後再加起這些乘積。以現性代數的話來說，「矩陣 M 用向量 amps 乘」。以訊號的話來說，我們計算了頻率成份的加權總和。
範例中，我標示第 n 個橫列用 a 到 d，它的值 會是 。
點積(dot product)的結果是 ys，它是一個向量每橫列是一個元素對應到 M。第 n 個元素標記為 e，是乘積的合：
ys 其他的元素也是如此計算。所以 y 的每個元素是四個頻率的總和，在某個時間點計算，且乘上對應的振幅。這正是我們要的。
我們可以用前一節的程式來驗證兩個版本的 synthesize 會得到一樣的結果。
ys1 = synthesize1(amps, fs, ts)
ys2 = synthesize2(amps, fs, ts)
max(abs(ys1 - ys2))寫下這個用線性代數的計算讓程式碼變小也變快。線性代數提供簡潔的記號給矩陣與向量。例如，我們可以將 synthesize 寫成這樣：
a 是振幅的向量，t 是時間的向量，f 是頻率的向量， 是兩個向量外積的記號。
6.3 分析 | Analysis
現在我們準備好來解決分析問題。假設我給你一個波，告訴你它是一組給定頻率的 cosine 的和。你該如何找到每個頻率成份的指幅？換句話說，給你ys、ts、fs，你能找回 amps 嗎？在線性代數的語言裡，第一步跟合成問題一樣：我們計算 。然後要找到一個 使得 。換句話說，我們要解出一個線性系統。NumPy 提供
linalg.solve，它正是做這事的。它的程式碼像這樣：
def analyze1(ys, fs, ts):
args = np.outer(ts, fs)
M = np.cos(PI2 * args)
amps = np.linalg.solve(M, ys)
return ampsnp.linalg.solve 計算 amps。這裡有個小技巧，一般我們只在 matrix 是正方形的時候解系統的線性方程組，也就是說方程式(rows)的數量要等同末知數(unknowns)的數量。
在這例子中，我們只有四個頻率，但我們評估這個訊號有 11025 次。所以我們的方程式遠比未知數多。
一般來說，如果
ys少於四個元素，我們無法只有四個頻率來分析它。但在這個例子，我們知道
ys只有四個頻率加總，所以我們可以用任意四個 wave array 的值來還原 amps。為了簡化，我會使用訊號的前面四個採樣，用前一節的
ys，fs，ts的值，來執行程式 analyze1：n = len(fs)
amps2 = analyze1(ys[:n], fs, ts[:n])amps2 會是 [0.6 0.25 0.1 0.05]這演算法有用，但是慢。解線性系統的時間會正比於 ， 是 的行數。我們可以做得更好。
6.4 正交矩陣 | Orthogonal matrices
一個解線性系統的方法是利用反矩陣。的反矩陣寫成 。它有個性質是 ， 是單位矩陣，它的對角線的元素全是 1，其他元素是 0。所以，要解 ，我們可以把兩邊都乘上 然後得到：
在右邊，我們把 換成 ：
如果 乘上向量 會得到 ，所以：
這暗示我們，如果能有效率地計算 ，我們就可以用簡單的矩陣乘法(用 np.dot)來找到 。這時間花費正比於 ，比 好。
一般找反矩陣是慢，但有些情況下比較快。一個特殊情況下，如果 是正交(orthogonal)， 的反矩陣就是 的轉置矩陣，寫成 。在 NumPy 中，轉置一個 array 是個常數時間的操作。它不是真的移動 array 的元素，取而代之的是建立一個視角，改變元素取用的方式。
一個矩陣要說是正交(orthogonal)，它的轉置矩陣要等同於反矩陣，也就是 ，這暗示 ，所以我們可以計算 來檢查是否一個矩陣是正交。
來看在 sysnthesize2 的矩陣長怎樣。在前一個例子裡， M 有 11025 列，所以用個小一點的例子會是個好主意。
def test1():
amps = np.array([0.6, 0.25, 0.1, 0.05])
N = 4.0
time_unit = 0.001
ts = np.arange(N) / N * time_unit
max_freq = N / time_unit / 2
fs = np.arange(N) / N * max_freq
ys = synthesize2(amps, fs, ts)amps我們之前看過，是振幅的向量表示。因為我們有四個頻率成份，我們會在訊號中的時間中取四個點。這樣 會是正方形。ts 也是等距時間間隔的時間取樣點的向量表示，範圍是 0 到 1 的時間單位。我選擇時間單位是 1ms。這是個可隨意的選擇，我們會看到選 1 分鐘算很慢就是了。因為 frame rate 是每時間單位 N 個採樣，Nyquist 頻率會是 ，這例子是 2000 Hz，所以
fs 是 0 到 2000 Hz 等距頻率的向量表示。有了這些
ts 與 fs 的值，矩陣 會是：[[ 1. 1. 1. 1. ]
[ 1. 0.707 0. -0.707]
[ 1. 0. -1. -0. ]
[ 1. -0.707 -0. 0.707]]但是，不是正交，如果我們計算 會得到：
[[ 4. 1. -0. 1.]
[ 1. 2. 1. -0.]
[-0. 1. 2. 1.]
[ 1. -0. 1. 2.]]6.5 DCT-IV
但是，如果我們小心選擇ts 與 fs，我們可以讓 成為正交。因為存在好幾個方法，所以也就有好幾個版本的 離散餘弦轉換 Discrete Cosine Transform (DCT)。一個最簡單的版本是移動
ts 與 fs 半個單位，這個版本叫做 DCT-IV，IV 是羅馬數字，代表它是 DCT 八個版本的第四個。現在來看改版過後的 test1：
def test2():
amps = np.array([0.6, 0.25, 0.1, 0.05])
N = 4.0
ts = (0.5 + np.arange(N)) / N
fs = (0.5 + np.arange(N)) / 2
ys = synthesize2(amps, fs, ts)這樣會得到的 M 是
[[ 0.981 0.831 0.556 0.195]
[ 0.831 -0.195 -0.981 -0.556]
[ 0.556 -0.981 0.195 0.831]
[ 0.195 -0.556 0.831 -0.981]][[ 2. 0. 0. 0.]
[ 0. 2. -0. 0.]
[ 0. -0. 2. -0.]
[ 0. 0. -0. 2.]]因為 是對稱，而且(幾乎)是正交， 的反矩陣就是 。現在我們可以寫更有效庠的
analyze：def analyze2(ys, fs, ts):
args = np.outer(ts, fs)
M = np.cos(PI2 * args)
amps = np.dot(M, ys) / 2
return ampsnp.linalg.solve，只改用乘上 。把
test2 與 analyze2 合併，我們可以寫出 DCT-IV 的實作：def dct_iv(ys):
N = len(ys)
ts = (0.5 + np.arange(N)) / N
fs = (0.5 + np.arange(N)) / 2
args = np.outer(ts, fs)
M = np.cos(PI2 * args)
amps = np.dot(M, ys) / 2
return amps如果我們搞定它，那這個函式可以解決分析問題。也就是，給定一個 ys，它應該可以回復 amps。我們可以如此測試：
amps = np.array([0.6, 0.25, 0.1, 0.05])
N = 4.0
ts = (0.5 + np.arange(N)) / N
fs = (0.5 + np.arange(N)) / 2
ys = synthesize2(amps, fs, ts)
amps2 = dct_iv(ys)
max(abs(amps - amps2))6.6 反離散餘弦轉換 | Inverse DCT
最後要注意，analyze2 與 synthesize2 幾乎相同，其中差異只有 analyze2 有把結果除以 2。我們可以用這觀察來計算「反離散餘弦轉換 | inverse DCT」：def inverse_dct_iv(amps):
return dct_iv(amps) * 2inverse_dct_iv 解決合成問題，它使用振幅的向量，然後回傳 wave array，也就是 ys。要測試它可用 amps 開始，使用 inverse_dct_iv 與 dct_iv，試試得到的與一開始的差別。amps = [0.6, 0.25, 0.1, 0.05]
ys = inverse_dct_iv(amps)
amps2 = dct_iv(ys)
max(abs(amps - amps2))6.7 DCT 類別 | The Dct class
圖6.2 400 Hz 的三角訊號，取樣率 10 kHz
thinkdsp 提供 Dct 類別，它封裝 DCT，就跟 Spectrum 類別封裝 FFT 一樣。要得到 Dct 物件，可以用 Wave 呼叫 make_dct。
signal = thinkdsp.TriangleSignal(freq=400)
wave = signal.make_wave(duration=1.0, framerate=10000)
dct = wave.make_dct()
dct.plot()make_dct 使用 DCT-II，這是最常用的 DCT 形式，由 scipy.fftpack 提供。
import scipy.fftpack
# class Wave:
def make_dct(self):
N = len(self.ys)
hs = scipy.fftpack.dct(self.ys, type=2)
fs = (0.5 + np.arange(N)) / 2
return Dct(hs, fs, self.framerate)Dct 提供 make_wave，它進行「反DCT」操作。我們可以這樣試試：
wave2 = dct.make_wave()
max(abs(wave.ys-wave2.ys))make_wave 使用 scipy.fftpack.idct：
# class Dct
def make_wave(self):
n = len(self.hs)
ys = scipy.fftpack.idct(self.hs, type=2) / 2 / n
return Wave(ys, framerate=self.framerate) 6.8 練習
For the following exercises, I provide some starter code in chap06starter.ipynb. Solutions are in chap06soln.ipynb.Exercise 1 In this chapter I claim that analyze1 takes time proportional to n3 and analyze2 takes time proportional to n2. To see if that’s true, run them on a range of input sizes and time them. In Jupyter, you can use the “magic command” %timeit.
If you plot run time versus input size on a log-log scale, you should get a straight line with slope 3 for analyze1 and slope 2 for analyze2.
You also might want to test dct_iv and scipy.fftpack.dct.
Exercise 2 One of the major applications of the DCT is compression for both sound and images. In its simplest form, DCT-based compression works like this:
- Break a long signal into segments.
- Compute the DCT of each segment.
- Identify frequency components with amplitudes so low they are inaudible, and remove them. Store only the frequencies and amplitudes that remain.
- To play back the signal, load the frequencies and amplitudes for each segment and apply the inverse DCT.
In order to make this method practical, you need some way to store a sparse array; that is, an array where most of the elements are zero. NumPy provides several implementations of sparse arrays, which you can read about at http://docs.scipy.org/doc/scipy/reference/sparse.html.
Exercise 3 In the repository for this book you will find a Jupyter notebook called phase.ipynb that explores the effect of phase on sound perception. Read through this notebook and run the examples. Choose another segment of sound and run the same experiments. Can you find any general relationships between the phase structure of a sound and how we perceive it?
沒有留言:
張貼留言