Octaveレッスン(11) - 自己相関から周波数を求める

*自己相関から周波数を求める
自己相関のグラフに見られる極値の周期は、元の信号の周期と合致する。
Octaveでは、findpeaks関数を使用して、極値のリストを得ることができる。

[pks, locs] = findpeaks(S,"DoubleSided", "MinPeakHeight", 0.3);

と実行すると、pksに±0.3以上の極値（極大値＆極小値）のリストが返され、locsに極値のインデクスが返される。
隣接する3つの極致のインデクスの差分を求めれば、周期が求まる。周波数は周期の逆数。

以下、極値の座標を求めて、周波数を求める関数の例
---- findCorrHz.m ----
function [WsHzVector, HzVector] = findCorrHz(R, Ws, Fs)
    if(nargin == 0)
      disp("usage: [WsHzVector, HzVector] = findCorrHz(R, Ws, Fs)");
      return
    end
    Thresh = 0.3;
    HzVector = [];
    WsHzVector = [];
    len=size(R)(1,1);
    for k = [1:Ws:len]
      if k+Ws > len
        break
      endif
      Sub=R(k:k+Ws-1);
      [pks, locs] = findpeaks(Sub,"DoubleSided", "MinPeakHeight", Thresh);
      locs_len = size(locs)(1,1);
      if locs_len <= 1
        Ts = -Ws;
        Hz = 0;
      else
        Ts = abs(locs(2,1) - locs(1,1));
        if locs_len == 2
          Ts += Ts;
        else
          Ts += abs(locs(3,1) - locs(2,1));
        endif
        Hz = Fs/Ts;
      endif
      HzVector=vertcat(HzVector, Hz);
      WsHzVector=vertcat(WsHzVector, repmat(Hz,1,Ws)');
    endfor     

endfunction
----------------

HzVectorには、自己相関をWsサンプルの区間ごとに求めた周波数が返る（長さR/Wsのベクトル）。

WsHzVectorには、求めた周波数を各サンプル区間ごとにWs個ずつ複製した、長さWsのベクトルが返る。

実行例：

[y, Fs] =audioread('Ongen2/egypt_charmera.wav');

Ws=1024, R=corrWs(y, Ws);

[WsHzVector, HzVector] = findCorrHz(R, Ws, Fs);

plot(R)

hold on, plot(y)

plot(WsHzVector/1000)

黄色：元の音声信号、青色：Wsサンプルごとに区切って求めた自己相関、紫色：周波数

参考：
  findpeaks 極大値の検出
  Octave - Forge - findpeaks