# 數位信號系統實習
###### tags: `matlab`

:::danger
>**因為PPT上傳了，所以有些東西有稍微改一下，也些打錯的地方都以刪除線刪除了。**
>[name=黃英睿][time=Sun, Jan 8, 2017 10:03 PM][color=#ff0000]
>
有發現錯誤和英睿說!!感謝!!
:::

[TOC]

# 1.雙變數回歸
題目：$y = w1x1+w2x2+b，已知x1,x2,y，求 w1,w2,b,在TwoVariables.mat裡面，各100,000筆。
## 程式碼

clear;
clc;
load('TwoVariables.mat');
%資料筆數
m = 100000;
%執行次數
time = 1000;
%預先宣告w1,w2,b
w1 = rand;
w2 = rand;
b = rand;
%設定步距
u = 0.1 ;
%設J的初始值
J = 0;
jw1 = 0;
jw2 = 0;
jb = 0;

for i = 1:time   
    J = 1/(2*m)*sum(((w1*x1+w2*x2+b)-y).^2)
    jw1 = mean(((w1*x1+w2*x2+b)-y).*x1);
    jw2 = mean(((w1*x1+w2*x2+b)-y).*x2);
    jb = mean((w1*x1+w2*x2+b)-y);
    w1 = w1 - (u*jw1);
    w2 = w2 - (u*jw2);
    b = b - (u*jb);
end

>*註解*
>sum()：加總
>mean()：取平均
## 結果(本題J的收斂值)
J =
2.0145

# 2.分析頻率響應並設計逆系統(並分析其頻率響應，理解方塊圖、差分方程式和Z轉換與辨識頻率響應之關係)
差分方程式 例: y[n]=(x[n]+x[n-1]+x[n-2])/3
Z轉換 例:
方塊圖 例:

將以上三樣求其頻率響應


H1為求得的頻率響應，H2則是其逆系統的頻率響應，兩者互為倒數

## 程式碼

clear;clc 
f=8000;
[y,fs]=audioread('test.mp4');
t=floor((size(y,1)/fs)*f);
tt=[1 round((1:t)/f*fs)];
for i=1:2
  yy(:,i)=y(tt,i);
end

d=2;
h=zeros(1,d*f);
h(1)=1;h(d*f)=0.5;
for i=1:2
  yn(:,i)=conv(yy(:,i),h);
end

Y=zeros(1,d*f);
Y(1)=1;Y(d*f)=0.5;
for i=1:2
  y1(:,i)=filter(1,Y,yn(:,i));
end

以聲音這題為例
這部分的差分方程式為 y[n]=x[n]+0.5x[n-D]
D為延遲的時間，即秒數d乘上取樣頻率f

d=2;
h=zeros(1,d*f);
h(1)=1;h(d*f)=0.5;
for i=1:2
  yn(:,i)=conv(yy(:,i),h);
end

而這部分的差分方程式則為 y[n]+0.5y[n-D]=x[n]

Y=zeros(1,d*f);
Y(1)=1;Y(d*f)=0.5;
for i=1:2
  y1(:,i)=filter(1,Y,yn(:,i));
end

---

f=8000;
[y,fs]=audioread('test.mp4');
t=floor((size(y,1)/fs)*f);
tt=[1 round((1:t)/f*fs)];
for i=1:2
  yy(:,i)=y(tt,i);
end

這一段是做把test.mp4的取樣頻率改成8k。

* 第 2 行：設定變數f(新的取樣頻率)。
* 第 3 行：做audioread這個方法，取出test.mp4中的y(聲音資料)和fs(test.mp4的取樣頻率)。
* 第 4 行：設定變數t，計算從fs(44.1K)換成f(8K)後共有多少點，取無條件捨去。
* 第 5 行：設定變數tt，把所有要從y中取出來值的位置都列出來，從1開始。(tt這個**陣列**後面的中括弧中放入兩個值，分別是一個**數字**，和一個**陣列**。大概是以下這種感覺)

a = [4 5 6];
b = [1 a];

    b = 
        [1 4 5 6]

* 第 6~8 行：做雙聲道的處理。
>註解:
>* [y,fs]=audioread(filename)：讀檔，y為讀出來的資料，fs為該資料的取樣頻率，filename為資料名稱。
>* audiowrite(y,fs,filename)：存檔，y為輸入來的資料，fs為該資料的取樣頻率，filename為資料名稱(***注意：matlab2013中如果要存成.mp4檔案的話，取樣頻率只能設44100~48000***)。
>* floor(a)：無條件捨去，將a做無條件捨去。
>* ceil(a)：無條件進位，將a做無條件進位。
>* round(a)：四捨五入，將a做四捨五入。
>[name=黃英睿][time=Mon, Jan 9, 2017 9:44 PM]
---

# 3.估計多路徑延遲
>題目：只給x[n]和y[n]求Delay。

解法：
1. 先猜有幾個Delay，先給Delay範圍(本題8000，P)，調整b~0~,b~1~,...,b~8000~ X[n-8000]。
2. $\mathrm{\hat{y}[n]-y[n] = e[n]}$，根據e[n]去做調整，其中$\mathrm{\hat{y}[n]}$是猜的。
3. 用filter把$\mathrm{\hat{y}[n]}$代入。
4. 計算方程式：
>$\mathrm{\hat{y}[n] = \sum^{P-1}_{k=0}b_{k}x[n-k]=b_{0}x[n]+b_{1}x[n-1]+...+b_{P-1}x[n-(p-1)]\\
>e[n]=y[n]-\hat{y}[n]=y[n]-\sum^{P-1}_{k=0}b_{k}x[n-k]\\
>E = (e[n])^{2}=(y[n]-\sum^{P-1}_{k=0}b_{k}x[n-k])^{2}\\
>\frac{dE}{dbm}=2×(y[n]-\sum^{P-1}_{k=0}b_{k}x[n-k])×(-x[n-m])=2e[n]x[n-m]\\
>~\\
>new\{b_{0},...,b_{P-1}\}\\
>~=old\{b_{0},...,b_{P-1}\}+u×2×e[n]×\{x[n],x[n-1],...,x[n-(P+1)]\}\\
>=>bm=bm+2×u(y[n]-\hat{y}[n])×x[n-m]}$

5. 把以上方程式轉換成以下程式碼

## Steepest Descend(梯度下降法)
是一個一階最佳化算法，通常也稱為最速下降法。 要使用梯度下降法找到一個函數的局部極小值，必須向函數上當前點對應梯度（或者是近似梯度）的反方向的規定步長距離點進行疊代搜索。如果相反地向梯度正方向疊代進行搜索，則會接近函數的局部極大值點；這個過程則被稱為梯度上升法。---[維基百科](https://zh.wikipedia.org/wiki/%E6%A2%AF%E5%BA%A6%E4%B8%8B%E9%99%8D%E6%B3%95)

## 程式碼

clear;
clc;
load('TwoDelay.mat');

P = 8000;
bm = zeros(1,P);
u = 0.001;

for xStart = 0:size(x,1)-P
    for loop = 1:100
        y_cap = bm*flipud(x(xStart+1:xStart+P));
        e = y(xStart+P)-y_cap;
        bm = bm + 2*u*e *flipud(x(xStart+1:xStart+P))';       
    end
    plot(bm)
    pause(0.00001)
end

11~13也可以用以下寫法，結果是一樣的。

        y_cap = filter(bm,1,x(xStart+1:xStart+P));
        e = sum(y(xStart+1:xStart+P)-y_cap);
        bm = bm + u * e.*flipud(x(xStart+1:xStart+P))';

>*註解*
>flipud()：將陣列反轉，例：A=[1 2 3]，flipud(A)=[3 2 1]。
>
>Y = filter(b,a,X)：卷積，X是卷積前的值x，Y是卷積後的值y，b是x的係數，a是y的係數。
>
>

#矩陣轉置
a = 
    [1 2 3]
a' = 
    1
    2
    3

## 程式解析

load('TwoDelay.mat');

* 第 3 行：load "TwoDelay.mat"裡面的參數，分別是x和y的值。
---

P = 8000;
bm = zeros(1,P);
u = 0.001;

* 第 5 行：設定參數的初始值P，P是猜的Delay的數量，猜8000是因為電磁波在真空中的傳輸速率大概是3×10^8^(m/s)，所以以取樣頻率8K來看，經過8000點後的值已經繞過了3×10^8^m了，這樣的值已經不值得參考了，所以取8000就已經足夠了。
* 第 6 行：設定參數的初始值bm，bm是我們要猜的x的係數，一開始都全部設為0。
* 第 7 行：設定參數的初始值u，u是步伐(step)，不能設太大，否則會無法收斂到最小值，但設太小就會跑很久。
---

for xStart = 0:size(x,1)-P
    for loop = 1:100

* 第 9 行：xStart是每次開始取x的位置，size是算x的陣列大小(本題是48112)，後面減P(本題是8000)是因為每次從xStart的位置向後取出P個點出來做運算，假設我們沒有減P這樣他會一直取到xStart = 48112(x的size)，然後向後取8000點，但是後面已經沒有點了，所以算出來會錯。對應第 17 行的end。
* 第 10 行：每次取P個點來做運算的次數。對應第 14 行的end。
---

        y_cap = bm*flipud(x(xStart+1:xStart+P));
        e = y(xStart+P)-y_cap;
        bm = bm + 2*u*e .*flipud(x(xStart+1:xStart+P))';

基本上就是把方程式寫成程式碼。
>方程式
>$\mathrm{\hat{y}[n] = b_{0}x[n]+b_{1}x[n-1]+...+b_{P-1}x[n-(p-1)]\\
>e[n]=y[n]-\hat{y}[n]\\
>b_m=b_m+2×u×e[n]×x[n-m]}$

* 第 11 行：這裡的 \* 不是一般的乘法，而是做矩陣的乘法。參考下面。
>bm =
>　　[b~0~ b~1~ ... b~P-2~ b~P-1~]
>x =
>　　x[1]
>　　x[2]
>　　$~~\vdots$
>　　x[n - 1]
>　　x[n]
>bm * x =
>　　b~0~x[1] + b~1~x[2] + ... + b~P-2~x[n-1] + b~P-1~x[n]
>bm * flipud(x) =
>　　b~0~x[n] + b~1~x[n-1] + ... + b~P-2~x[2] + b~P-1~x[1]

* 第 12 行：照方程式的打。
* 第 13 行：照方程式的打，最後面 分號(**;**) 前有一個 單引號(**'**)， 單引號(**'**)可以把這矩陣做矩陣轉置。做這行的目的是為了讓bm後面那串(u * e.\*flipud(x(xStart+1:xStart+P)))的陣列型態和bm相同，這樣才可以做矩陣的加法。

---

    end
    plot(bm)
    pause(0.00001)
end

* 第 14 行：對應第 10 行的for
* 第 15~16 行：把每次做完就把bm plot出來看結果，pause是暫停時間，讓程式有時間把圖畫出來。
* 第 17 行：對應第 9 行的for

## 結果
約跑6分鐘後，在大約1200及4800有明顯的值。

# 4.偵測頻率orDTMF訊號(內含混音雜訊)

## 偵測頻率
### 程式碼

clear;
clc;

f1 = 3e3;
f2 = 5e3;
fs = 20e3;
t = 0:1/fs:0.1;
N = 2^10
f=(1:N)/N*fs;

x = cos(2*pi*f1*t)+cos(2*pi*f2*t);
y = fft(x,N);
plot(f,abs(y));
axis([0 10000 200 550]);

>*註解*
>fft(x,N)：快速傅立葉轉換，可以取出此x的頻率，N是取樣點，通常設成2的指數，但不能超過x的點數(**出來的值是複數**)。
>abs()：取絕對值，在複數上是取大小，例 ：a = 3+4j，abs(a)=5；b = -1，abs(b) = 1。
>第 9 行：f=(1:N)/N\*fs，是讓1到N的刻度變成1到fs的刻度。
>這題如果沒有打第 14 行也可以，會看到4個地方有值，如下圖，但後面兩個是不對的(老師說的但是忘記是為什麼了XD)，只要看前面兩根的頻率就好了。
>

### 結果

## DTMF
>**雙音多頻信號**（**D**ual-**T**one **M**ulti-**F**requency, **DTMF**），電話系統中電話機與交換機之間的一種用戶信令，通常用於發送被叫號碼。
>
>在使用雙音多頻信號之前，電話系統中使用一連串的斷續脈衝來傳送被叫號碼，稱為脈衝撥號。脈衝撥號需要電信局中的操作員手工完成長途接續。---[維基百科](https://zh.wikipedia.org/wiki/%E5%8F%8C%E9%9F%B3%E5%A4%9A%E9%A2%91 "雙音多頻 - 維基百科，自由的百科全書")

### DTMF.mat中Tone內容

題目：~~總時間10秒~~，取樣頻率(Fs)=8kHZ，共有10個號碼，每個號碼持續0.5秒並間隔0.5秒。
>刪除線處說明:這裡應該是不會給總時間，但可以自己算出來，總時間的算法是總點數÷取樣頻率(Fs)[name=黃英睿][time=Sun, Jan 8, 2017 10:00 PM]

解法：
1. 找出總點數：~~10(總時間(秒))\*8k(每秒取8k點)＝80k(80,000)點~~ 載入DTMF.mat後看裡面Tone的大小。

2. 每個號碼持續0.5秒並間隔0.5秒，所以每個號碼有 0.5(持續時間)\*8K(取樣頻率) = 4000個點，每次只要取4000點出來做FFT分析出頻率。

### 程式碼

clear;
clc;
load('DTMF.mat');

fs = 8e3;
N = 2^11;
for a =0:9
    f = (1:N)/N*fs;
    T=Tone(fs*a+1:fs*a+4001);
    y = fft(T,N);
    ny = abs(y);
    subplot(5,2,a+1);
    plot(f,ny);
    axis([500 2000 0 1500]);
end

>*註解*
>axis([Xmin Xmax Ymin Ymax])，設定plot中xy軸的顯示範圍。

### 結果
對照下表得出：0277387411
### DTMF表
|\\|1209HZ|1336HZ|1447HZ|1633Hz|
|:--:|:--:|:--:|:--:|:--:|:--:|
|**697HZ**|1|2|3|A|
|**770HZ**|4|5|6|B|
|**852HZ**|7|8|9|C|
|**941HZ**|\*|0|#|D|

# 5.設計濾波器(低通、高通、帶通、帶斥)，包含脈衝響應(h)與頻率響應(H)
**題目：** 以20K作為取樣頻率(Fs)設計出以下四個Filter
1. 截止頻率(Fc)為**3k**的**Low-Pass Filter**
2. 截止頻率(Fc)為**5k**的**High-Pass Filter**
3. 截止頻率(Fc)範圍為**3k~5k**的**Band-Pass Filter**
4. 截止頻率(Fc)範圍為**3k~5k**的**Band-Rejection Filter**

## Filter
* Low-Pass Filter：$h_{lp}[n] = \frac{\Omega_{c}}{\pi}sinc(\frac{\Omega_{c}}{\pi}(n-L))$
* High-Pass Filter：$h_{hp}[n] = \delta[n-L] - \frac{\Omega_{c}}{\pi}sinc(\frac{\Omega_{c}}{\pi}(n-L))$
* Band-Pass Filter：$h_{bp}[n] = \frac{\Omega_{c1}}{\pi}sinc(\frac{\Omega_{c1}}{\pi}(n-L)) *(\delta[n-L] - \frac{\Omega_{c2}}{\pi}sinc(\frac{\Omega_{c2}}{\pi}(n-L)))$
* Band-Rejection Filter：$h_{br}[n] = \delta[n-L] - (\frac{\Omega_{c1}}{\pi}sinc(\frac{\Omega_{c1}}{\pi}(n-L)) *(\delta[n-L] - \frac{\Omega_{c2}}{\pi}sinc(\frac{\Omega_{c2}}{\pi}(n-L)))$

## Other Windowing Effect
>在信號處理中，Window function（窗口函數）是一個數學函數是零值的一些選擇的外面間隔。例如，一個功能是時間間隔內常數和零別處被稱為矩形窗口，其描述了它的圖形表示的形狀。當另一個函數或波形/數據序列乘以窗口函數時，乘積也在間隔之外為零值：剩下的是它們重疊的部分，即“通過窗口查看”。
>
>在典型的應用中，所使用的窗口函數是非負的平滑的“鐘形”曲線。矩形，也可以使用三角形，以及其它功能。窗口函數的更一般的定義，不要求它們是相同的零的時間間隔以外，只要在窗口乘以其參數的乘積為平方可積，並且，更具體地，該功能是為足夠迅速趨向零。---[維基百科](https://en.wikipedia.org/wiki/Window_function "Window function - Wikipedia")
>

Hanning：$w(k) = 0.5 - 0.5cos(2\pi\frac{k}{N})$，for k=0,...,N-1
Hamming：$w(k) = 0.54 - 0.46cos(2\pi\frac{k}{N})$，for k=0,...,N-1
Blackman：$w(k) = 0.42 - 0.5cos(2\pi\frac{k-1}{N})+0.08cos(4\pi\frac{k-1}{N})$，for k=0,...,N-1

### 程式碼

clear;
clc;

Fs=20e3;

OmegaC1 = pi/((Fs/2)/3e3); %3K
OmegaC2 = pi/2; %5K
L = 20;
n = 0:2*L;

N1 = 2*L+1;
k1 = 0:N1-1;
N2 = 4*L+1;
k2 = 0:N2-1;
Hanning1 = 0.5 - 0.5*cos(2*pi*k1/N1);
Hanning2 = 0.5 - 0.5*cos(2*pi*k2/N2);

dat = zeros(1,2*L+1);
dat(L+1) = 1;
dat2 = zeros(1,4*L+1);
dat2(2*L+1) = 1;

h_lp = (OmegaC1/pi*sinc(OmegaC1/pi*(n-L)));
h_hp = dat - (OmegaC2/pi*sinc(OmegaC2/pi*(n-L)));
h_bp = conv((OmegaC2/pi*sinc(OmegaC2/pi*(n-L))),dat - (OmegaC1/pi*sinc(OmegaC1/pi*(n-L))));
h_br = dat2 - h_bp;

[Hlp,w]=freqz(h_lp .* Hanning1);
wf = w/pi*Fs/2;
subplot(2,2,1)
plot(wf,abs(Hlp));
title('低通，Ωc=3K');

[Hhp,w]=freqz(h_hp .* Hanning1);
wf = w/pi*Fs/2;
subplot(2,2,2)
plot(wf,abs(Hhp));
title('高通，Ωc=5K');

[Hbp,w]=freqz(h_bp .* Hanning2);
wf = w/pi*Fs/2;
subplot(2,2,3)
plot(wf,abs(Hbp));
title('帶通，Ωc1=3K、Ωc2=5K');

[Hbr,w]=freqz(h_br .* Hanning2);
wf = w/pi*Fs/2;
subplot(2,2,4)
plot(wf,abs(Hbr));
title('帶斥，Ωc1=3K、Ωc2=5K');

>*註解*
>~~截止頻率(fh)的算法(假設截止頻率是5k、取樣頻率為20K)： $OmegaC = \frac{\pi}{\frac{f_{sh}}{f_{h}}} = \frac{\pi}{\frac{10kHZ}{5kHZ}} = \frac{\pi}{2}$~~
>上面的方程式打錯了，雖然答案會一樣，但是以下才是正確的寫法(來源:老師的PPT)[name=黃英睿][time=Sun, Jan 8, 2017 8:45 PM]
>$\Omega$c的算法：fc(截止頻率，假設5k)、fs(取樣頻率，假設20k)
>$OmegaC = fc\frac{2\pi}{fs} = 5k\frac{2\pi}{20k} = \frac{\pi}{2}$
>

OmegaC2 = pi/2;

>顯示頻率，取樣頻率的一半：$f_{sh} = \frac{fs}{2}$
>>[H,w]=freqz(h)，計算**h**(脈衝響應)的**H**(頻率響應，為複數)和**W**(角頻率，值為0~$\pi$)[color=#ff0000]

### 結果脈衝響應圖(H)(沒有使用Windowing Effect)
### 結果脈衝響應圖(H)(使用Hanning、Hamming和Blackman)