電子科技大學(xué)數(shù)字信號處理DSP課程設(shè)計-鋼琴音符識別
.
2014級數(shù)字信號處理課程設(shè)計報告
題目：鋼琴音符識別
姓名：邱晨曦
學(xué)號：2014010909008
答辯時間：2016/12/9
一. 題目要求：
（1） 播放和記錄一段鋼琴音樂中的音符；
（2） 記錄到音符以后，找到音符所對應(yīng)的現(xiàn)代標(biāo)準(zhǔn)鋼琴的鋼琴鍵，并分析結(jié)果。
二. 課程設(shè)計思路：
（1） 涉及到的知識點：
快速傅里葉變換、鋼琴音頻信號的時域和頻域的特性、能熵比的概念、頻率校正、頻率與音符的轉(zhuǎn)換關(guān)系。
（2） 方案分析：
A. 預(yù)處理部分：
1. 直接用audioread函數(shù)讀出來的原始數(shù)據(jù)。
優(yōu)點：準(zhǔn)確率較高；
缺點：數(shù)據(jù)量較大，采樣頻率為44kHz，遠(yuǎn)大于奈奎斯特采樣率。
2. 以11kHz的采樣率重新采樣，并轉(zhuǎn)換為單聲道。
優(yōu)點：數(shù)據(jù)量小了很多，易于處理；
缺點：犧牲了部分的準(zhǔn)確率，但對于音符的判斷影響可以忽略。
B. 端點檢測算法：
<1> . 雙門限法：
1. 計算短時能量（高門限）和過零率（低門限）；
2. 選取一個較高的門限，語音信號的能量包絡(luò)大部分都在此門限之上，進行一次初判，語音起止點位于該門限與短時能量包絡(luò)交點所對應(yīng)的時間間隔之外；
3. 根據(jù)噪聲能量，確定一個較低的門限，并從初判起點往左，從初判終點往右搜索，分別找到能零比曲線第一次與門限相交的兩個點，兩點之間段就是用雙門限方法所判定的語音段；
4. 以短時平均過零率為準(zhǔn)，從低門限點往左右搜索，找到短時平均過零率低于某閾值的兩點，為語音的起止點；
圖1：雙門限法示意圖
說明：算法中的閥值是根據(jù)實驗過程調(diào)節(jié)的。
該算法在實際應(yīng)用的過程中發(fā)現(xiàn)：在語音信號頻率分布較為集中的時候，端點檢測出來的結(jié)果比較準(zhǔn)確，但當(dāng)語音信號頻率分布比較分散的時候，很難通過控制固定的閥值來檢測到每個音符；
<2>. 自相關(guān)法：
由于兩種信號的自相關(guān)函數(shù)存在極大的差異，可以利用這種差別來提取語音端點。根據(jù)噪聲的情況，設(shè)置兩個閾值和，當(dāng)相關(guān)函數(shù)最大值大于時，便判定是語音；當(dāng)相關(guān)函數(shù)最大值大于或小于時，則判定為語音信號的端點。
該算法同樣存在當(dāng)語音信號頻率分布較廣的時候，閥值比較難控制的問題。
<3>. 基于譜熵的端點檢測：
基于譜熵語音端點檢測方法是通過檢測譜的平坦程度，來進行語音端點檢測的，為了更好進行語音端點檢測，采用語音信號的短時功率譜構(gòu)造語音信息譜熵，從而對語音段和噪聲段進行區(qū)分。檢測思路：
1. 對語音信號進行分幀加窗；
2. 計算每一幀的譜能量；
3. 計算出每一幀中每個樣本點的概率密度函數(shù)；
4. 計算出每一幀的譜熵值（由信息論知識知道，熵值在自變量服從均勻分布的時候，熵值達(dá)到最大值，所以噪聲的熵值是比較大的，而鋼琴音符的熵值是比較小的，由此區(qū)別了噪聲和音符）；
5. 設(shè)置判決門限；
6. 根據(jù)各幀的譜熵值進行端點檢測。
在實驗過程中發(fā)現(xiàn)：依然存在當(dāng)語音信號頻率分布較廣時，閥值不太好控制的問題。因此對該方法進行改進，引入，能熵比的概念：
譜熵值類似于過零率，能熵比的表示為。由于噪聲和信號的能熵比差別很大。因此在能熵比的圖像中，每一個“尖刺”就代表了一個特定頻率的語言信號。
圖2：能熵比圖中的“尖刺”
在檢測過程中，依然不能通過簡單的設(shè)置閥值的辦法來進行端點檢測，原因是語音頻率分布較廣時，每個音符的能熵比變化范圍差別較大，如下圖所示，有的“尖刺”完全在門限之上，而有的則完全在門限之下。
圖3：88階全音的能熵比圖
因此，采用檢測能熵比中的“低谷點”（該點比左右兩邊的一定數(shù)目的點的能熵比都?。┑姆椒?。語音信號一定位于兩個低谷點之間的部分，再對低谷點進行適當(dāng)?shù)淖笥乙苿幼鳛檎Z音信號的起止點。如下圖所示：
圖4：標(biāo)記起止點的能熵比圖
（綠色為起始點，紅色為截止點）
（3） 設(shè)計框架和流程：
1. 用audioread函數(shù)讀入鋼琴音樂，并用sound函數(shù)播放；
2. 為了方便處理，對信號以11.025kHz的頻率進行重新采樣，并統(tǒng)一轉(zhuǎn)換成單聲道的信號；
3. 因為語言信號可以在短時間內(nèi)認(rèn)為是平穩(wěn)的，因此對語音信號進行分幀的處理，設(shè)置幀長320，為了減小誤差，兩幀之間設(shè)置重疊部分，因此幀移取80；
4. 計算每一幀的能熵比；
5. 找到能熵比中的“低谷點”（該點比左右兩邊的一定數(shù)目的點的能熵比都小）；
6. 如果兩個低谷點之間的距離大于miniL(認(rèn)為持續(xù)長度超過一定長度的為音符，最小長度miniL可自行設(shè)置)。則低谷點右移sr (即shift right,數(shù)值可自行調(diào)節(jié))幀作為一段信號的起始點，將低谷點左移sl（即shift left,數(shù)值可自行調(diào)節(jié)）幀作為截止點[注：采用該方法的優(yōu)點是通過調(diào)節(jié)相關(guān)參數(shù)能適應(yīng)多種情況，缺點是檢測環(huán)境發(fā)生較大變化時，需要重新設(shè)置參數(shù)]；
7. 將找到的語音段轉(zhuǎn)換成未分幀時對應(yīng)坐標(biāo)的語音段，并對每段做快速傅里葉變換；
8. 找到每段快速傅里葉變換中的最大值以及最大值所對應(yīng)的橫坐標(biāo)（fft點），將橫坐標(biāo)轉(zhuǎn)換成相應(yīng)的頻率，得到的頻率即為該段音符的頻率；
9. 利用比值法進行頻率的校正，窗函數(shù)選擇矩形窗；
10. 根據(jù)檢測到的頻率確定音符，計算公式為：，為第幾個按鍵，再通過查表得到對應(yīng)音符；
11. 分析結(jié)果。
三. 具體設(shè)計過程：
（1） 部分代碼（測試部分缺?。?
主函數(shù)部分：
[x,fs]=audioread(鋼琴音頻.WAV);
format short;
wlen=320; inc=80; % 分幀的幀長和幀移
overlap=wlen-inc; % 幀之間的重疊部分
sound(x,fs); % 播放音樂
x=calsample(x,fs); % 為了方便處理，重新以11025Hz的頻率采樣，并轉(zhuǎn)換成單聲道
x=x-mean(x); % 消去直流分量
x=x/max(abs(x)); % 幅值歸一化
y = Enframe(wlen,inc,x); % 分幀
fn = size(y,2); % 取得幀數(shù)
time = (0 : length(x)-1)/11025; % 計算時間坐標(biāo)
frameTime = frame2time(fn, wlen, inc, 11025); % 計算各幀對應(yīng)的時間坐標(biāo)
sr=2;sl=13;miniL=33; % 配置左右移動的幀數(shù)和要求的最短幀數(shù)
[voicesegment,vos,Ef]=get_segment(y,fn,sr,sl,miniL); % 獲得語音段
[real_f,ft,ax]=get_f(x,voicesegment,vos,wlen,inc); % 檢測頻率的結(jié)果
for i=1:length(real_f)
real_f(i)=roundn(real_f(i),-4);
end
for i=1:length(real_f)
real_node(i)=get_node(real_f(i))
end
%**********************************繪圖部分********************************
subplot 211;
stem(real_f);
title(頻率檢測結(jié)果);xlabel(音符/個);ylabel(頻率/Hz);
subplot 212;
stem(real_node,r);
title(音符檢測結(jié)果);xlabel(音符/個);ylabel(對應(yīng)按鍵);
figure(2);
subplot 211;
plot(Ef);
title(能熵比圖及語音起止點);xlabel(幀數(shù)/個);ylabel(能熵比);
for i=1:length(voicesegment)
text(voicesegment(i).begin,Ef(voicesegment(i).begin),o,color,g)
text(voicesegment(i).end,Ef(voicesegment(i).end),o,color,r)
end
subplot 212, plot(time,x,k); title(語音信號端點檢測結(jié)果)
axis([0 max(time) -1 1]); ylabel(幅值);
for k=1 : vos % 標(biāo)出有話段
nx1=voicesegment(k).begin;
nx2=voicesegment(k).end;
nxl=voicesegment(k).duration;
fprintf(%4d %4d %4d %4d\n,k,nx1,nx2,nxl);
subplot 212
line([frameTime(nx1) frameTime(nx1)],[-1 1],color,r,linestyle,-);
line([frameTime(nx2) frameTime(nx2)],[-1 1],color,r,linestyle,--);
end
其中的用到的子函數(shù)：
1.calsample.m （調(diào)整采樣率和聲道）
function sample = calsample(sampledata,FS)
temp_sample = resample(sampledata,1,FS/11025); %調(diào)整采樣頻率
[~,n] = size(temp_sample);
if (n == 2) %轉(zhuǎn)換成單聲道
sample = temp_sample(:,1);
else
sample = temp_sample;
end
end
2. Enframe.m （分幀函數(shù)）
function f=Enframe(len,inc,x) %對讀入的語音進行分幀，len為幀長，
%inc為幀重疊樣點數(shù)，x為輸入語音數(shù)據(jù)
fh=fix(((size(x,1)-len)/inc)+1); %計算幀數(shù)
f=zeros(fh,len); %設(shè)置一個零矩陣，行為幀數(shù)，列為幀長
i=1;
n=1;
while i<=fh %幀間循環(huán)
j=1;
while j<=len %幀內(nèi)循環(huán)
f(i,j)=x(n);
j=j+1;
n=n+1;
end
n=n-len+inc; %下一幀開始位置
i=i+1;
end
3. frame2time.m（坐標(biāo)刻度轉(zhuǎn)換）
function frameTime=frame2time(frameNum,framelen,inc,fs)
frameTime=(((1:frameNum)-1)*inc+framelen/2)/fs; % 求對應(yīng)的時間坐標(biāo)
4.get_segment.m（端點檢測，確定音符段）
function [voicesegment,vos,Ef]=get_segment(y,fn,sr,sl,miniL)
if size(y,2)~=fn, y=y; end % 把y轉(zhuǎn)換為每列數(shù)據(jù)表示一幀語音信號
wlen=size(y,1); % 取得幀長
for i=1:fn
Sp = abs(fft(y(:,i))); % FFT取幅值
Sp = Sp(1:wlen/2+1); % 只取正頻率部分
Esum(i) = sum(Sp.*Sp); % 計算能量值
prob = Sp/(sum(Sp)); % 計算概率
H(i) = -sum(prob.*log(prob+eps)); % 求譜熵值
end
hindex=find(H<0.1);
H(hindex)=max(H);
Ef=sqrt(1 + abs(Esum./H)); % 計算能熵比
Ef=Ef/max(Ef); % 歸一化
[x1,y1] = get_max(Ef); % 找到能熵比中的“高峰點”
[x2,y2] = get_min(Ef); % 找到能熵比中的“低估點”
voicesegment(1).begin = x1(1)-8; % 由于僅僅靠低谷點無法檢測出第一個音符起始位置， 因此將高峰點的第一個值左移8幀作為第一個音符的起始點
voicesegment(1).end = x2(1)-sl; % 將第一個低谷點作為第一個音符的截止點
voicesegment(1).duration=voicesegment(1).end -voicesegment(1).begin+1; % 一個音符的持續(xù)幀數(shù)
j=1;
for k=2 : length(x2) % 將找到的低谷點作為音符的起止點
if x2(k)-x2(k-1)>=miniL % 剔除持續(xù)幀長度小于miniL的音符段
j=j+1;
temp1=x2(k-1)+sr; % 將低谷點右移sr個幀，作為一個音符的起始點
voicesegment(j).begin=temp1;
temp2=x2(k)-sl; % 將低谷點左移sl個幀，作為一個音符的截止點
voicesegment(j).end=temp2;
voicesegment(j).duration=voicesegment(j).end-voicesegment(j).begin+1; %音符持續(xù)幀數(shù)
end
end
vos=length(voicesegment); %返回音符個數(shù)
end
5.get_max.m（找到高峰點）
function [max_x,max]=get_max(x)
l=length(x); %獲得數(shù)組的長度
max=[];
max_x=[];
j=1;
for i=20:l-15 %找到“峰值點”
if (x(i)>x(i-1))&&(x(i)>x(i-2))&&(x(i)>x(i-3))&&(x(i)>x(i-4))&&(x(i)>x(i-5))...
&&(x(i)>x(i-6))&&(x(i)>x(i-7))&&(x(i)>x(i-8))&&(x(i)>x(i-9))&&(x(i)>x(i-10))...
&&(x(i)>x(i-11))&&(x(i)>x(i-12))&&(x(i)>x(i-13))&&(x(i)>x(i-14))...
&&(x(i)>x(i-15))&&(x(i)>x(i-16))&&(x(i)>x(i-17))...
&&(x(i)>x(i+1))&&(x(i)>x(i+2))&&(x(i)>x(i+3))&&(x(i)>x(i+4))&&(x(i)>x(i+5))...
&&(x(i)>x(i+6))&&(x(i)>x(i+7))&&(x(i)>x(i+8))&&(x(i)>x(i+9))&&(x(i)>x(i+10))...
&&(x(i)>x(i+11))&&(x(i)>x(i+11))&&(x(i)>x(i+12))&&(x(i)>x(i+13))&&(x(i)>x(i+14))...
&&(x(i)>x(i+15))&&(x(i)>x(i+16))...
&&(x(i)>0.1)
max_x(j)=i; %找到后賦值給返回參數(shù)
max(j)=x(i);
i=i+1;
j=j+1;
else
i=i+1;
end
end
end
6.get_min.m（找到低谷點）
function [min_x,min]=get_min(x)
l=length(x); %獲得數(shù)組的長度
min=[];
min_x=[];
j=1;
for i=100:l-10 %尋找低谷點
if (x(i)<x(i-1))&&(x(i)<x(i-2))&&(x(i)<x(i-3))&&(x(i)<x(i-4))&&(x(i)<x(i-5))...
&&(x(i)<x(i-6))&&(x(i)<x(i-7))&&(x(i)<x(i-8))&&(x(i)<x(i-9))&&(x(i)<x(i-10))...
&&(x(i)<x(i-11))&&(x(i)<x(i-12))&&(x(i)<x(i-13))&&(x(i)<x(i-14))...
&&(x(i)<x(i-15))&&(x(i)<x(i-16))&&(x(i)<x(i-17))&&(x(i)<x(i-18))&&(x(i)<x(i-19))...
&&(x(i)<x(i-20))&&(x(i)<x(i-21))&&(x(i)<x(i-22))&&(x(i)<x(i-23))&&(x(i)<x(i-24))...
&&(x(i)<x(i-25))&&(x(i)<x(i-26))&&(x(i)<x(i-27))&&(x(i)<x(i-28))&&(x(i)<x(i-29))...
&&(x(i)<x(i-30))&&(x(i)<x(i-31))&&(x(i)<x(i-32))&&(x(i)<x(i-33))&&(x(i)<x(i-34))......
&&(x(i)<x(i+1))&&(x(i)<x(i+2))&&(x(i)<x(i+3))&&(x(i)<x(i+4))&&(x(i)<x(i+5))...
&&(x(i)<x(i+6))&&(x(i)<x(i+7))&&(x(i)<x(i+8))&&(x(i)<x(i+9))&&(x(i)<x(i+10))...
min_x(j)=i; %找到后賦值給返回值
min(j)=x(i);
i=i+1;
j=j+1;
else
i=i+1;
end
end
end
7.get_f.m（確定頻率）
function [real_f,ft,ax]=get_f(x,voicesegment,vos,wlen,inc)
for i=1:vos % 橫坐標(biāo)轉(zhuǎn)化，將起止點坐標(biāo)的單位由幀轉(zhuǎn)換成點
ax(i).begin=(voicesegment(i).begin-1)*inc+1;
ax(i).end=(voicesegment(i).end-1)*inc+wlen;
ax(i).duration=ax(i).end-ax(i).begin+1;
temp=x(ax(i).begin:ax(i).end); %獲得語音段
ft{i}=fft(temp); % 做快速傅里葉變換，并將結(jié)果保存在元胞矩陣中
ft{i}=abs(ft{i}); % 將得到的fft去模值
[fm(i),fm_x(i)]=max(ft{i}); % 找到每個語音段對應(yīng)的fft中的最大值及最大值對應(yīng)的橫坐標(biāo)
real_f1(i)=fm_x(i)*11025./ax(i).duration; % 進行頻率轉(zhuǎn)換，將橫坐標(biāo)乘以分辨率得到真實頻率
real_f(i)=Specorrm(x(ax(i).begin:ax(i).end),11025,ax(i).duration,real_f1(i)-10,real_f1(i)+10);
% 利用比值法進行頻率修正，窗函數(shù)選擇矩形窗
end
8.Specorrm.m（頻率校正）
function Z=Specorrm(x,fs,N,nx1,nx2)
%x是被測信號, fs是采樣頻率 N為FFT的長度，nx1和nx2被測信號頻率的區(qū)間，nx2>nx1
[nx,mx]=size(x);
if mx==1, x=x;end %轉(zhuǎn)換成行矩陣
[nx,mx]=size(x);
M=fix(N/2)+mod(N,2);
xf=fft(x);
%xf=xf(1:M)*2/N;
ddf=fs/N; % 頻率分辨率
n1=fix(nx1/ddf); % 將頻率轉(zhuǎn)換成fft對應(yīng)的點
n2=round(nx2/ddf);
A=abs(xf); % 取fft的模值
[Amax,index]=max(A(n1:n2)); %找到fft點n1到n2之間幅值的最大值
index=index+n1-1; % 移動到n1和n2中間
%比值法
%加矩形窗
indsecL=A(index-1)>A(index+1); %滿足條件則為1，不滿足則為0
df=indsecL.*A(index-1)./(Amax+A(index-1))-(1-indsecL).*A(index+1)./(Amax+A(index+1));
Z=(index-1-df)*ddf; %修正后的頻率
end
9.get_node.m（確定音符）
function node=get_node(f)
node=12*log2(f/440)+49;
node=round(node); %四舍五入確定對應(yīng)按鍵
end
（2） 生成的圖與數(shù)據(jù)等
端點檢測結(jié)果
序號
起始幀
截止幀
持續(xù)幀數(shù)
1
40
97
58
2
112
166
55
3
181
235
55
4
250
303
54
5
318
372
55
6
387
441
55
7
456
510
55
8
525
579
55
9
594
648
55
10
663
714
52
11
729
786
58
12
801
854
54
13
869
923
55
14
938
993
56
15
1008
1051
44
16
1066
1122
57
17
1137
1176
40
18
1191
1267
77
19
1282
1337
56
20
1352
1391
40
21
1406
1475
70
22
1490
1513
24
23
1528
1588
61
24
1603
1666
64
25
1681
1730
50
26
1760
1808
49
27
1823
1873
51
28
1903
1944
42
29
1959
2017
59
30
2032
2090
59
31
2105
2136
32
32
2177
2225
49
33
2240
2278
39
34
2317
2340
24
35
2380
2429
50
對數(shù)據(jù)的分析：
真實按鍵
真實音符
真實頻率/Hz
檢測頻率/Hz
對應(yīng)按鍵
檢測音符
絕對誤差/Hz
相對誤差
52
G?7/A?7
523.251
524.8973
52
G?7/A?7
1.6463
0.31%
53
G7
554.365
555.0437
53
G7
0.6787
0.12%
54
F?7/G?7
587.33
589.8768
54
F?7/G?7
2.5468
0.43%
55
F7
622.254
623.1585
55
F7
0.9045
0.15%
56
E7
659.255
659.7079
56
E7
0.4529
0.07%
57
D?7/E?7
698.456
699.2968
57
D?7/E?7
0.8408
0.12%
58
D7
739.989
741.8503
58
D7
1.8613
0.25%
59
C?7/D?7
783.991
785.9485
59
C?7/D?7
1.9575
0.25%
60
C7Double high C
830.609
832.414
60
C7Double high C
1.805
0.22%
61
B6
880
883.2303
61
B6
3.2303
0.37%
62
A?6/B?6
932.328
933.8565
62
A?6/B?6
1.5285
0.16%
63
A6
987.767
989.9964
63
A6
2.2294
0.23%
64
G?6/A?6
1046.5
1047.01
64
G?6/A?6
0.51
0.05%
65
G6
1108.73
1112.7244
65
G6
3.9944
0.36%
66
F?6/G?6
1174.66
1177.1137
66
F?6/G?6
2.4537
0.21%
67
F6
1244.51
1248.1238
67
F6
3.6138
0.29%
68
E6
1318.51
1322.609
68
E6
4.099
0.31%
69
D?6/E?6
1396.91
1403.0084
69
D?6/E?6
6.0984
0.44%
70
D6
1479.98
1486.3162
70
D6
6.3362
0.43%
71
C?6/D?6
1567.98
1571.9648
71
C?6/D?6
3.9848
0.25%
72
C6Soprano C(High C)
1661.22
1669.9962
72
C6Soprano C(High C)
8.7762
0.53%
73
B5
1760
1773.532
73
B5
13.532
0.77%
74
A?5/B?5
1864.66
1881.0102
74
A?5/B?5
16.3502
0.88%
75
A5
1975.53
1988.0017
75
A5
12.4717
0.63%
76
G?5/A?5
2093
2112.8655
76
G?5/A?5
19.8655
0.95%
77
G5
2217.46
2241.0368
77
G5
23.5768
1.06%
78
F?5/G?5
2349.32
2377.199
78
F?5/G?5
27.879
1.19%
79
F5
2489.02
2512.4893
79
F5
23.4693
0.94%
80
E5
2637.02
2515.4901
79
F5
121.5299
4.61%
81
D?5/E?5
2793.83
2824.2907
81
D?5/E?5
30.4607
1.09%
82
D5
2959.96
2994.1387
82
D5
34.1787
1.15%
83
C?5/D?5
3135.96
3192.2427
83
C?5/D?5
56.2827
1.79%
84
C5Tenor C
3322.44
3364.7355
84
C5Tenor C
42.2955
1.27%
85
A7
3520
3575.9524
85
A7
55.9524
1.59%
86
A?7/B?7
3729.31
3771.6582
86
A?7/B?7
42.3482
1.14%
紅色標(biāo)記出檢測錯誤的音符。
上表為檢測結(jié)果的統(tǒng)計及誤差統(tǒng)計。大部分誤差都小于1%，平均誤差為0.82%，誤差較小.音符檢測發(fā)生一個錯誤，正確了為97.14%，很好地完成了任務(wù)。
四. 總結(jié)：
（1） 收獲
（2） 小組分工情況
小組就邱晨曦一個人，邱晨曦負(fù)責(zé)全部工作。
參考文獻(xiàn)：謝明，丁康：離散頻譜的一種新校正方法，重慶，載重慶大學(xué)學(xué)報：1995年3月 第8卷第2期 47~54
.
個人主頁