clear,clc
% 1. 對給定語料估計其基音周期。
% 要求用MATLAB或C語言實現(xiàn)有關(guān)基音檢測算法，并給出檢測結(jié)果。
% 
% 2. 實驗語料：MaoYiSheng.raw
% 發(fā)音：“茅以升”；時長0.9s；采樣率：8kHz；量化精度為16bits/sample；數(shù)據(jù)格式 LSB,MSB
% 分析幀長:20ms(160樣點)
% 總幀數(shù)：45   45*160=7200
% 3. 兩周內(nèi)將算法描述、源程序、結(jié)果及其分析通過E-MAIL發(fā)送到 wangjing@bit.edu.cn
N=160;
fid = fopen('MaoYiSheng.raw', 'r');
% Y1= fread(fid);
Y = fread(fid,'short');
figure(1),
% subplot(4,1,1),
plot(Y);
title('原始譜');
xlabel('采樣點個數(shù)');

grid;
L=length(Y);
LL=length(Y)/N;
% energy
Em=zeros(1,(LL-1)*160);
for ii=1:(LL-1)*160
    temp=Y(ii:ii+160-1);
    temp=temp.*hamming(160);
    Em(ii)=sum(temp.*temp);
end
figure(2),
% subplot(4,1,2)
jj=[1:(LL-1)*160];
plot(jj,Em);
title('短時能量');
xlabel('采樣點個數(shù)');
ylabel('E');
grid,
length(Y);
for j=1:89
    %  j=50;   
    y=Y(((j-1)*80+1):(j+1)*80);
    x=y.*boxcar(160);
    % R=zeros(1,160);
    % for k=1:160
    %     for n=1:160-k
    %         R(k)=R(k)+x(n)*x(n+k);
    %     end
    % end
    r=zeros(1,160);
    for k=1:160
        for n=1:160-k
            r(k)=r(k)+x(n)*x(n+k);
        end
    end
    % figure(3),
    % subplot(4,1,3);
    % j=1:160;
    % plot((j-1)*80:(j+1)*80-1,R);
    % title('濁音自相關(guān)');
    % xlabel('滯后點個數(shù)k');
    % ylabel('相關(guān)R');
    % grid;
    r=r(21:160);
    [Z,I]=max(r);
    R(j)=I+21;
end
figure(3),
% subplot(4,1,3);
plot(1:80:89*80,R);
title('周期圖');
xlabel('語音采樣點');
ylabel('周期長度');
grid;