Octaveレッスン(12) - 区間ごとに異なる周波数の正弦波を生成する

* 区間ごとに異なる周波数の正弦波を生成する

「Octaveレッスン(11) - 自己相関から周波数を求める」で区間ごとに周波数（HzVector）を求めた。

以下、HzVectorで与えられる周波数情報を元に、区間ごと正弦波を生成する例
---- genWsSin.m ----

function WsSin = genWsSin(HzVector, Ws, Fs)
    if(nargin == 0)
      disp("usage: WsSin = genWsSin(HzVector, Ws, Fs)");
      return
    endif
    len=size(HzVector)(1,1);
    preHz=-1;
    dur=0;
    WsSin=[];
    for k = [1:len-1]
      crrHz=HzVector(k,1);
      if crrHz==preHz
        dur += Ws;
      else
        if dur > 0
          WsSin = horzcat(WsSin, gen_pre_sin(preHz, dur, Fs));
        endif
        preHz=crrHz;
        dur=Ws;
      endif      
    endfor

    if dur > 0
      WsSin = horzcat(WsSin, gen_pre_sin(preHz, dur, Fs));
    endif
endfunction

function y = gen_pre_sin(preHz, dur, Fs)
  x=[0:dur-1];
  y=sin(2*pi*preHz/Fs*x);

endfunction

----------------

実行例：

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

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

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

WsSin = genWsSin(HzVector, Ws, Fs);

sound(WsSin,Fs)

⇒ 笛の音を録音した音声信号を1024サンプルの区間に分割して、区間ごとに周波数を検出。

　 genWsSin()関数を使用して、区間ごとに正弦波で再構成した音WinSinを再生すると、元の音色に近い音階が聞こえる

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

Octaveレッスン(10) - 自己相関を求める

*自己相関を求める

音声の周波数を求めるには、自己相関が用いられることが多い。
自己相関のグラフに現れる極大値の周期が、音声の周波数に該当する。

Octaveでは、xcorr() 関数を使用して、自己相関を簡単に求めることができる。

以下、音声信号SをWsサンプルずつに分割して、Wsサンプルごとに自己相関を求める関数の例

---- corrWs.m ----

function r = corrWs(S,Ws)
  if(nargin == 0)
    disp("usage: r = corrWs(S, Ws)");
    return
  end

  %% if stereo, use only 1ch.
  if size(S)(1,2) > 1
    S = S(:,1);
  endif
  len = size(S)(1,1);
  r = [];
  for k = [1:Ws:len]
    if (k+Ws-1) > len
      return
    end
    s1=S(k:k+Ws-1);
    r1 = xcorr(s1, s1, 'coeff');
    r = vertcat(r, r1(1:Ws-1));
  end
endfunction

----------------

【NOTE】
・xcorr() の引数に 'coeff' オプションを指定すると正規化した値が返される。xcorr() に 'coeff' 引数を指定して自己相関を求める場合には、同じ関数同士の相互相関として、引数を指定する必要がある。
・自己相関は左右対称なので、前半のみ使用している

実行例：

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

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

plot(R)

hold on, plot(y)

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

参考： xcorr 相互相関

Octaveレッスン(9) - FIRフィルタについて復習

*FIRフィルタについて復習

昔、FIRフィルタについて、お勉強したハズなのだが、すっかり忘れてしまったので、参考書籍をもろもろ読み直して復習した。

ここに要点を記載しておく（また忘れないように）。

まずは、f(t)のフーリエ変換、逆フーリエ変換の式をおさらい：

  F(ω)=∫f(t)e^(-jωt) dt　　…(1)
  f(t)=1/2π∫F(ω)e^(jωt)dt　　…(2)

上記2式は下式のように略記されることが多い：

  F(ω)=F{f(t)}
  f(t)=F^-1{F(ω)}

線形システムにおいて、インパルス応答h(t)のフーリエ変換をH(ω)とすると、下式が成り立つ：

　f(t)＊h(t)=1/2π∫F(ω)H(ω)e^(jωt)dt

時間領域の畳み込み演算（＊）だが、周波数領域では単なる積の演算になっている。

なお、線形システムとは、入力信号が足し算で表されるなら出力信号もその足し算で表され、入力が定数倍になれば、出力も定数倍になるような系のこと。

ローパスフィルタ（周波数ωc以下の低域を通し、それより上の周波数は遮断）は、下図の左上のような方形窓H(ω)を周波数領域でF(ω)と乗算することを意味する。

方形窓H(ω)を逆フーリエ変換したh(t)を求めて、時間領域にて、h(t)と入力信号f(t)の畳み込み演算を行うと、ローパスフィルタのかかった出力信号が得られることになる。

方形窓の逆フーリエ変換を計算するとsinc関数になる：

 h(t)=1/πt sin(ωt)

実用の場では、離散フーリエ変換を使用するので、下図の右上のような離散的なsinc関数になる。

sinc関数は無限に続く関数なので、このままでは使えない（有限時間で計算できない）。
そのため、1/πt sin(ωn) のnをどこかで打ち切ってフィルタ係数として使用する（下図の中央右）。

『はじめて学ぶディジナル・フィルタと高速フーリエ変換』（CQ出版）p.93 図6-1

有限範囲で係数を打ち切ると、上図の中央左のように、通過域・阻止域にうねり（リップル）が生じてしまう（ギブス現象）。
フィルタ係数の打ち切りは、時間領域で矩形窓をかけていることを意味する。
リップルを和らげるために、時間領域でなだらかな窓関数（ハミング窓など）と係数を乗算する（上図の右下）。

得られたフィルタ係数の周波数特性を求め（上図の左下）、所望の特性が得られていることを確認する。

以上でFIRフィルタについての復習はおしまい。

ついでにフーリエ変換の対称性についても復習しておく。

(1), (2)を見比べると、2式には類似性がある。
(2)において、tを-ω, ωをtと置き換えると、

 2πf(-ω)=∫F(t)e^(-jωt)dt
              = F{F(t)}            

つまり、フーリエ変換の対称性を示す下式が得られた：

  F{F(t)} = 2πf(-ω)

フーリエ変換の対称性が役立つのは、f(t)のフーリエ変換F(ω)が判明しているが、F(ω)のωをtに入れ替えた関数F(t)のフーリエ変換が簡単に求まらない場合。

F(t)のフーリエ変換は、対称性からf(t)のtを-ωに入れ替えたf(-ω)に2πを乗じた関数になる。

例えば、方形窓r(t)

r(t) = 1(-1<t<1), 0(t<-1 or 1<t)

のフーリエ変換R(ω)を求めると

R(ω) =∫1*e^(-jωt)dt = 2 sin(ω)/ω

これを踏まえると、関数f(t):

f(t) = sin(ω)/ω

のフーリエ変換を求めることができる。
f(t)=R(t)/2なので、フーリエ変換の対称性より、

F(ω) = 2πr(-ω)/2 = πr(-ω)

と簡単に求まる。

Octaveレッスン(8) - バンドパスフィルタをかける

*バンドパスフィルタをかける

butter関数を使用して各種バタワースフィルタを作成できる。
標準の関数ではないので、事前にsignalパッケージをロードするコマンドの実行が必要：

pkg load signal

次数10、低域カットオフ周波数 440 Hz 、高域カットオフ周波数 660 Hzのバンドパスフィルターを作成する例：

Fs = 44100;     % sampling rate
[b, a] = butter(10, [440 660]/(Fs/2));

[b, a] に、伝達関数の係数が返される。

500Hzの正弦波と800Hzの正弦波の合成波に、バンドパスフィルタをかける：

[y500, t] = gen_sin(Fs, 500); % sampling rate 44100Hz, sin 500Hz
[y800, t] = gen_sin(Fs, 800); % sampling rate 44100Hz, sin 800Hz
yy = y500 + y800;
sound(y500, Fs)  % play sound
sound(y800, Fs)  % play sound
sound(yy, Fs)  % play sound
yy2=filter(b, a, yy);
sound(yy2, Fs)  % play sound

合成波yyから800Hzの成分が除去されて、yy2には500Hzの成分だけ残る。

sound()コマンドでyy2を再生するとy500と似た再生音が聞こえる。

・・・となるはずだが、yy2が破裂音のような変な音になってしまう。

先日、同様の手順でうまくいったのだが要調査。

参考： バタワース フィルターの設計

Octaveレッスン(7) - AndroidでOctaveを使う

* Android版Octaveのセットアップ方法

Google Playから次の3つのアプリをインストールするだけでよい。
3.をインストールすることで、plotコマンドも使えるようになる。
plot画像をpngファイルとして生成してくれる。

1. Octave
2. Octave Main Package
3. Octave Gnuplot Package

毎起動時にアップデートをお勧めされるが、私の旧めの環境（Android4.2.2）では、新バージョンのインストールに失敗した。
現行バージョンで、軽快に何の問題もなく動作してくれるので "Maybe Later" を選択。

* Android版Octaveが使用するディレクトリ

私の端末で使用されていたディレクトリ（環境によって異なる可能性がある）：

• scriptファイル(.m, .octaverc)置き場： /storage/emulated/0/GNUOctave/home
• plotファイル（.png）置き場： /storage/emulated/0/GNUOctave/intents

試しに、Octaveレッスン(4) - 正弦波を生成する関数 の.mファイルを作成して、

[y, t] = gen_sin(Fs, hertz, seconds);

実行すると、wavファイル作成には失敗したが（audiowrite関数がない）、出力 y, t を得られた。
plot(t,y) コマンドも実行できた。
ただし、Android端末の画面サイズが小さいので、hertzが小さくないと、波形が潰れてしまい、青い矩形が描画されただけのpng画像が生成される。