信号処理のお仕事メモ

2025年6月22日日曜日

Python GUIアプリ入門講座「世界にひとつだけのPython GUIアプリ作成入門」キャンペーン情報 ★2025/7/22 まで！★

久しぶりですが、Udemyキャンペーン情報をお知らせします

世界にひとつだけのPython GUIアプリ作成入門

※ 上記リンクからキャンペーン価格 1600円で購入可能（2025/7/22 まで）

クーポンコード： 6BB09DF6FCC0FA48B8C7

2022年7月31日日曜日

もっとも簡単なPythonインストール方法とGUIアプリ作成方法

Pythonインストールするもっとも簡単な方法は、Microsoft Storeを使う方法ではないでしょうか。

あまり知られていない気がするので、動画でご紹介：

そして、GUIアプリを作成するもっとも簡単な方法は、PySimpleGUIを使う方法だと思うのですが、これまたなぜか、あまり知られていません。

公式ページ： PySimpleGUI

主に英語ですが（いちおう日本語のreadmeもあります）、サンプル満載で詳しく解説されています。

有償ですがUdemy講座だとGUIの動きを確認しながら動画で学べるので、分かりやすいかも。

開発者Mike氏の公式Udemy講座： Python GUIs - "The Official PySimpleGUI Course"
※ 上記リンクからキャンペーン価格 1730円で購入可能（2022/12/29 現在）

日本語で学びたい方向け：世界にひとつだけのPython GUIアプリ作成入門
※ 上記リンクからキャンペーン価格 1600円で購入可能（2023/1/29 まで）

PySimpleGUIは、凝ったデザインのGUIアプリを作るには不向きですが、少ないコード量で、簡単にGUIアプリを作成できます。

コンソールアプリだと、ユーザーに敬遠されてしまったり、使用感が悪くなってしまうということがありがちです。PySimpleGUIを使って、ちょっとだけ手間をかけてGUIアプリ化してあげると、アプリの可能性が大きく広がります。

2021年5月24日月曜日

シェーダーで遊ぼう！

シェーダーは、GPUで実行される、描画の並列処理を記述するためのプログラムです。

CPUで実行するためのGDScriptとは別物。
GDScriptはPythonベースですが、シェーダーはC言語ベースの異なる言語です。

シェーダーの種類(shader_type)には、spacial(3D用), canvas_item(2D用), particles(パーティクルシステム用)の3種類があります。まずは canvas_item を使いこなせるようになりましょう。

CPUで実行する、横width, 縦heightの矩形領域を赤色に設定する処理をC言語で書くと、以下のようなコードになります。
for文で、矩形領域内の座標(x, y)のピクセルを順番に処理。
set_col() 関数で、座標(x, y)のピクセルの色を red に設定するというイメージです。

for(int x = 0; x < width; x++) {
for(int y = 0; y < height; y++) {
set_col(x, y, red);
}
}

シェーダーの世界では、上記のようなfor文が無くなります。
並列に実行される、描画領域内の各ピクセルにおける処理のみを記述するので、以下のようなコードになります。

// function called for each pixel
void fragment() {
  COLOR = red;
}

COLORは、ピクセルの色を表す変数。
redは、R,G,B, アルファで構成される４次元ベクトルで、
red = vec4(1.0, 0.0, 0.0, 1.0);
という値です。

fragment() 関数は、プロセッサー関数と呼ばれるもので、他に vertex() 関数、light() 関数という二人の仲間がいます。一番よく使うのは、色の描画処理を行う fragment() 関数なので、まずは、この関数と仲良くなりましょう。

シェーダーを使うには、ColorRectなどのnodeのインスペクタで、Materialに「新規 Shader Material」を選択。
shaderが「空（empty）」となっているところに、「新規 Shader」を選択すると、コード編集画面が開きます。

英語版ですが、以下の操作動画がわかりやすいです（kids can codeより引用）。

シェーダーのコードを書いてみましょう。

shader_type canvas_item;

void fragment() {

COLOR = vec4(UV.x, UV.y, 0.5, 1.0);

}

これを実行すると、以下のようにColorRectノードの色が変わります。

(UV.x, UV.y)は「描画領域の左上端を(0, 0)、右下端を(1, 1)」とみなした座標系における座標で、正規化座標と呼ばれるものです。正規化座標を使うことで、実際の描画領域のサイズが異なる場合でも、同じコードを利用できます。

上記のシェーダーコードでは、RGBのRをUV.xにしているので、右に行くほど赤くなります。
GはUV.yなので、下に行くほど緑になります。Bは0.5の固定値なので、どこでも一定です。
vec4(UV.x, UV.y, 0.5, 1.0) の最後要素の1.0は、アルファ値という透明度を表すものです。
1.0は不透明、0.0は透明になります。

私が作成したペンライトアプリには、ライトを点滅させるモードがあります。
これを実現しているのは、次のシェーダーコードです。

shader_type canvas_item;

void fragment() {

vec3 red = vec3(1.0, 0.0, 0.0);

float gain = abs(sin(5.*TIME));

if(gain > 0.4) gain = 1.0;

COLOR = vec4(red * gain, 1.0);

}

TIMEは経過時間ごとに与えられる秒数です。
赤色の度合いをsin関数の絶対値を取って、0.0 〜1.0 の範囲で、時間変化させています。
光っている時間は、しっかり赤色が見えるようにしたいので、0.4以上の時は赤色の度合いを 1.0 に切り上げるようにしています。

もう一つのグラデーションモードも、30行ほどのシェーダーコードで実現しています。
皆さんも、自分なりのペンライトの光らせ方をさせてみると、よい練習問題になるとかもしれません。

ぜひシェーダーで遊んでみてください！

楽しいですよ〜

ペンライトぴかっち

■Godot入門記事

1. Godot EngineでAndroid/iOSアプリを同時開発しよう！
2. Godot Engineおすすめ学習コース
3. Godot Engineの入門チュートリアルが終わったら
4. シェーダーで遊ぼう！

2021年5月13日木曜日

Godot Engineの入門チュートリアルが終わったら

Godot Engineのオフィシャル入門チュートリアルステップ・バイ・ステップが終わったら、いよいよオリジナルゲームを作ってみましょう！

といきたいところですが、作る準備は整いましたか？
私はダメでした。。。

手順に従って作成したゲーム「Dodge the Creeps」を自分のスマホで動かすことができたので、もう怖いものなし。どんなものでも作れそうな気分です。
しかし、いざ作ってみようとすると、手が止まってしまいます。。。

いきなり凝ったオリジナルアプリを作ろうとするのは、まだ無理でした。
そこで、ハードルを下げて、単純だけど使えるアプリを何か作れないかと考えて、ペンライトアプリを作ってみることにました。

コンサート会場でお馴染みのペンライトをスマホで実現しちゃおう！というアプリです。
画面遷移なし、セーブデータなしという超シンプル設計ですが、いちおう実用には耐えうるかと。

ペンライトぴかっち

5/9(日)にリリース申請して、やっと審査を通過。
本日5/13(木)より公開中です！

作り方については、のちほど別の機会にご紹介したいと思いますが、何でもよいので超シンプルなアプリを自分で考えて作ってみることをお勧めします。
もやっと理解していたことが、きちんと理解出来るようになったり、理解が足りないところが浮き彫りになります。

それと並行して、他のチュートリアルにも取り組んでみました。
まだまだ分かっていないことも多いですが、少しずつ知識の幅が広がっていると感じます。

私のお勧めは Godot Recipes by KidsCanCodeですが、他にも色々な学習サイトがあります。

オフィシャルサイトでも紹介されています：
チュートリアルとリソース

■Godot入門記事

2021年5月12日水曜日

Godot Engineおすすめ学習コース

Godot の入門チュートリアルは、日本語で丁寧に書かれているのですが、そのまま忠実になぞって学習を進めると、終了までにかなりの時間を要します。

まずは手っ取り早く動かして、どのようなアプリやゲームが作れるのか、知りたいところ。
さまざまなサンプルアプリが用意されていますので、ダウンロードして遊んでみましょう！

Godot オフィシャルページAsset Library で Category に project を指定すると、色々なサンプルプロジェクトを見つけることができます。

こちらにも各種デモプロジェクトが取り揃えられています： godot-demo-projects

■実行手順

Godotエディタを起動して、「インポート」ボタンを押して、ダウンロード＆展開したしたフォルダ内のproject.godt というファイルを指定します。

インポートしたプロジェクトを選択してダブルクリック（または「実行」ボタン押し）。
Godotエディタの編集画面が開きます。

サンプルを実行するには、右上の三角ボタンを押すか、
ショートカットキー：
・windowsなら F5
・Macなら command+B
を入力。

次は、サンプルを少しだけカスタマイズして、遊んでみましょう。

いきなりゲームをちょい変するのは、ハードルが高いので、
godot-demo-projects の「GUI」フォルダ内にある、GUIパーツのサンプルプロジェクトが手頃でお勧め。

例えば、Drag＆Drop というサンプルがあります。

色の付いた四角形をドラッグして、他の四角形へドロップすると、投げ込まれた色に変化するというアプリです。これをちょい変して、投げ込まれた色と元の色をブレンドした色に変化するようにしてみましょう。

シーンタブの内の「ColorPickerButton0」などの右横に、フニャっと斜めになった白い四角いアイコンがあるので、そこにカーソルを当てると「スクリプトを開く： ...」と表示されます。

そのままクリックすると、ボタンに割り当てられているスクリプトが開いてくれます。

あるいは、スクリプトのファイル名は「*.gd」という決まりがあるので、「ファイルシステム」タブ内の drag_drop_script.gd を見つけて、ダブルクリックしてもO.K.

細かい話は置いておいて、コード（たったの18行！）を眺めてみると

func drop_data(_pos, data):
color = data

のところで色を変更しているようなので、ここのコードを書き換えてみましょう。

元の色(color)と投げ込まれた色(data)を2で割るように、

func drop_data(_pos, data):
color = (color + data) / 2.0

と書き換えます。

実行してみてください。
思いどおりの動きに変更できましたね！

どんなものが作れるのか、ようすが分かってきたら、本腰を入れて、

オフィシャル入門チュートリアル
ステップ・バイ・ステップ

に取り組みましょう。時間はかかりますが、やる価値「大」です。

チュートリアルの最初のゲームまで終えると、ちょっとしたゲームが完成。
自分のスマホで動かすこともできます！
ディアゴスティーニで最終巻まで買い続けて、模型が完成したときのような感動が得られますヨ！！

ただ、結構な長丁場になりますので、ときどき気分転換することは重要かと。

・作りたい自作アプリ・ゲームの案を考えてみる（妄想していると、やる気が出ます！）
・またサンプルプロジェクトのちょい変をしてみる
・もし可能なら、ちょっとしたオリジナルアプリを作ってみる（無理は禁物）

といった気分転換を織り交ぜつつ、気長に取り組んでみてください！

ちなみに、私の場合、時間の取れるときに、少しずつ入門チュートリアルを進めていったのですが、２ヶ月ぐらいかかりました。
時間が空くと忘れてしまい、手戻りが発生しがちなので、ちょっとだけでもよいので毎日継続するよう心がけました。

■Godot入門記事

2021年5月10日月曜日

Godot EngineでAndroid/iOSアプリを同時開発しよう！

Godot Engineをご存知でしょうか？

あまり知名度は高くないのですが、無償で使えるゲームエンジンです。

ゲームエンジンというとUnityやUnreal Engineが有名です。しかし、多機能過ぎて、習得に時間がかかったり、ハイスペックなPCが必要になるという嫌いがあります。
実のところ、どちらにも挑戦したことがあるのですが、挫折しました。

そんなときGodot Engineを見つけました。
機能が多過ぎず、見通しがよいと感じました。あれこれデザインパターンで縛りつけようとせず、自由度の高い開発ができそうな雰囲気も気に入りました。
ファイルサイズも小さく、ダウンロードしたexeを叩けばそのまま起動するという手軽さも魅力。

WindowsやMacで作ったアプリを Android / iOS 用に変換してくれる機能もあります。

3D, 2Dゲームいずれも作成可能。2Dゲーム開発についてのみ比較すると、UnityやUnreal Engineより、機能が充実しているのではないでしょうか。

密かに、音声信号処理を利用した2Dゲーム的なモバイルアプリを作ってみたいという野望をもっている私は、今年の1月からGodot Engineの習得に取り組むことにしました。

またもや挫折しそうになったときもありましたが、なんとか持ち堪えて、このGWに自作アプリ第一号が完成！
ただいま、Google Playリリースに向けて審査中です。

第一号アプリは、「まずは何かオリジナルアプリを作ってリリースする」が目標で、かなり単純なアプリではあるのですが、とにもかくにも自分のアイデアを形にすることができました！嬉しい！！

今後の記事で、アプリ完成までの道のりを書いて行きたいと思います。
Godot Engineおすすめですヨ！

■ダウンロードはこちらから： Godot download
　※ C#信者以外は、standard を選択

■日本語チュートリアルも充実：ステップ・バイ・ステップ

■Godot入門記事

2017年1月1日日曜日

Octaveレッスン(16) - surf関数で3次元影付き表面プロット

明けましておめでとうございます。
surf関数で、新年の3Dグラフ表示をしてみました。

下記 bird_3d 関数を実行すると、ニワトリ風な物体が表示されます。
surf関数で複数の楕円体を表示しています。

---- bird_3d.m ----

function bird_3d
t = 0:0.05*pi:pi;
p = 0:0.05*pi:2*pi;
figure
%body
a=b=c=10
x = a*sin(t)'*cos(p);
y = b*sin(t)'*sin(p);
z = c*cos(t)'*ones(size(p));
surf(x,y,z)
xlabel('x'),ylabel('y'),zlabel('z')
hold on
%eyes
a=b=c=1
x = a*sin(t)'*cos(p);
y = b*sin(t)'*sin(p);
z = c*cos(t)'*ones(size(p));
surf(x+9,y+3,z+3)
surf(x+9,y-3,z+3)
%foot
a=8,b=2,c=1
x = a*sin(t)'*cos(p);
y = b*sin(t)'*sin(p);
z = c*cos(t)'*ones(size(p));
alpha=pi/8;
x_org=x; y_org=y;
y=cos(alpha)*y_org - sin(alpha)*x_org;
x=sin(alpha)*y_org + cos(alpha)*x_org;
surf(x+2,-y+4,z-9)
surf(x+2, y-4,z-9)
%crown
a=b=2,c=4;
x = a*sin(t)'*cos(p);
y = b*sin(t)'*sin(p);
z = c*cos(t)'*ones(size(p));
surf(x+2,y,z+9)
surf(x,y,z+10)
surf(x-2,y,z+9)
%bill
a=6,b=2,c=1
x = a*sin(t)'*cos(p);
y = b*sin(t)'*sin(p);
z = c*cos(t)'*ones(size(p));
surf(x+10,y,z)
surf(x+10,y,z-1)
%wings
t = 0:0.1*pi:pi;
p = 0:0.1*pi:2*pi;
a=4,b=8,c=1
x = a*sin(t)'*cos(p);
y = b*sin(t)'*sin(p);
z = c*cos(t)'*ones(size(p));
alpha=-pi/8;
y_org=y;z_org=z;
y=cos(alpha)*y_org - sin(alpha)*z_org;
z=sin(alpha)*y_org + cos(alpha)*z_org;
alpha=-pi/8;
x_org=x;z_org=z;
x=cos(alpha)*x_org - sin(alpha)*z_org;
z=sin(alpha)*x_org + cos(alpha)*z_org;
surf(x, y+8,z-1)
surf(x,-y-8,z-1)
endfunction

----

参考：

2016年12月26日月曜日

Octaveレッスン(15) - plot3関数で3次元グラフ表示

plot3関数を使用すると３次元グラフを簡単に表示することができる。

plot3(x, y, z) のようにx、y、zそれぞれの座標データを引数に与えればよい。

Octaveレッスン(6) で見た、FFTの基底の原始N乗根 ωN を３次元グラフに表示してみる。

ここでは N=128とする。

>> N=128, omega=e^(i*(2*pi/N));
>> x=0:N*20-1;
>> f=omega .^ x;
>> plot3(x, real(f), imag(f))
>> grid on
>> xlabel('X'), ylabel('real'), zlabel('imag')

グラフウィンドウの回転アイコン（左上端の円状の矢印）をクリックすると、グラフを回転させて、さまざまな角度から眺めることができる。
ぐるぐるグラフを回転させて x-y(real)平面を表示させると、cos関数のグラフが現れる。
さらに回転させて x-z(imag)平面を表示させると、sin関数のグラフが現れる。
y(real)-z(image)平面では、単位円のグラフが現れる。

参考：

甲南大学 Matlab/Octave関連情報第5章プロット応用
MathWorks 線形 3 次元プロット plot3

2016年12月24日土曜日

Octaveレッスン(14) - ltfat packageを利用する

ltfat packageには、fwt関数などwavelet解析に便利な関数が用意されている。

windows版Octaveでは、下記コマンドで ltfat のインストールに成功するが、 Mac版では失敗する。

>> pkg install -forge ltfat

ドキュメントを読む限り、Mac版もWindows版も同じ手順でインストールできないとおかしいのだが、別のMacや旧いOSバージョン、旧いOctaveバージョンでも同様に失敗する。

そこで、別のアプローチを試してみたところ成功したので、メモを残しておく。

Mac版で動かすには、下記ページからltfat packageをダウンロードして展開する：
http://ltfat.sourceforge.net/download.php

Octaveのコマンドウィンドウにて、展開した ltfat ディレクトリへ移動して、

>> ltfatstart

と実行すれば、各種 ltfat の関数が利用可能になる。
waveletのデモを表示するには、

>> demo_wavelets

さらに処理速度を高速化したい場合には、

>> ltfatmex

と実行すればよい・・・ハズなのだが失敗する。
fftw3ライブラリが見つからず、コンパイルエラーが発生する。

下記URLからfftw-3.3.5.tar.gzをダウンロードしてfftw3ライブラリをコンパイルする。
http://www.fftw.org/download.html

ターミナルを開いて、展開したfftw-3.3.5ディレクトリに移動して下記コマンドを実行：

$ ./configure --enable-threads

$ make

libfftw3.a, libfftw3_threads.a が生成されるので、場所をfind等で探して、ltfat/lib へコピーする。

続いて、float対応版ライブラリをコンパイルする：

$ make clean

$ ./configure --enable-float --enable-threads

$ make

libfftw3f.a, libfftw3f_threads.a が生成されるので、場所をfind等で探して、ltfat/lib へコピーする。

Octaveのコマンドウィンドウに戻って、

>> ltfatmex

と再実行すると、今度は成功する。

※ ltfatmex 実行すると、PortAudio が無いというコンパイル警告が出るため、私はPortAudioもダウンロード＆インストールした。しかし、これは必須では無いはず。

　PartAudio関連サイト：

　　http://www.portaudio.com/download.html

　　http://portaudio.com/docs/v19-doxydocs/compile_mac_coreaudio.html

Octaveレッスン(13) - signal packageを利用する

fft関数などは、Octaveをインストールするだけで利用可能。
しかし、xcorr関数（相互相関）などを利用しようとすると、signal packageが必要になる。

windows版の場合は、Octaveインストール時にsignalパッケージもインストールされるので、下記コマンドでsignal packageをロードすれば利用可能になる

>> pkg load signal

Macの場合は、まずinstallする必要がある：

>> pkg install -forge signal

上記は、ネットワーク経由でダウンロード＆インストールするコマンドだが、おそらくcontrol package が未インストールのため失敗するだろう。

そこで、まず control packageをinstallする：

>> pkg install -forge control

依存関係のあるcontrolをインストールできたので、signalをインストールしてロードする：

>> pkg install -forge signal

>> pkg load signal

これで、xcorr関数などが利用可能になる。

なお、インストール済みパッケージは下記コマンドで一覧表示可能：

>> pkg list

また、別 package 無しで利用可能な関数の一覧は、[ドキュメント]タブ - 「Function Index」sectionで確認可能。

2016年11月27日日曜日

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を再生すると、元の音色に近い音階が聞こえる

2016年11月26日土曜日

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 相互相関

2016年11月15日火曜日

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(-ω)

と簡単に求まる。

2016年11月5日土曜日

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ファイルとして生成してくれる。

Octave
Octave Main Package
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画像が生成される。

2016年10月31日月曜日

Octaveレッスン(6) - FFTのイメージを掴む

ωN を１の原始N乗根とすると、離散フーリエ変換の式は、以下のように書ける。

『これなら分かる応用数学教室』 p.130

離散フーリエ変換は、サンプル数Nのサンプルデータf0, ..., fN-1をN次元のベクトルとみなし、基本正弦波の0倍～N-1倍の高調波で直交展開することを意味している。

K≠Lのとき、K倍の高調波とL倍の高調波は常に直交する（p.114）。

FFT（高速フーリエ変換）とは、離散フーリエ変換のNを2のべき乗にすることで、フーリエ係数の計算を効率化する方式。
そうすることで、フーリエ係数を求める際に、ωN^kが計算過程で重複して登場するようになり、いちど算出した結果を使いまわすことで、計算量を減らすことができる。

サンプルfが実数のとき、フーリエ係数F0, ..., FN-1のうち、後ろ半分については、前半を折り返した複素共役になっていて、独立ではない(p.116)。
　∵)ωN^(N-k)l = conj(ωN^kl )なので、fが実数のときF(N - k) ｌ＝conj(Fｋl)

そのため、これまでOctaveで振幅スペクトルを求めてプロットする際は、プロットする範囲を0..N/2に半減し、F0, ..., FN/2 の絶対値を求めて2倍していた（省略した後半の分を加算するので2倍）。

フーリエスペクトル（フーリエ係数）： Fk
振幅スペクトル： Fkの絶対値 |Fk|
位相スペクトル： Fkの偏角 ∠Fk
エネルギースペクトル： |Fk|^2

参考書籍：『これなら分かる応用数学教室』p.129～p.139

2016年10月30日日曜日

Octaveレッスン(5) - 数式を描く

* 数式を描く

離散フーリエ変換の式をLatex形式で記述し、描画してみた。

Latexの全ての記述に対応している訳ではない。

分数を記述する\frac, \cfrac は使えなかった。

>> clf
>> text(0.1, 1-0.1, 'f_{i}=1/N \Sigma_{i=0}^{N-1} F_{k}\omega_{N}^{kl}')
>> text(0.1, 1-0.2, 'F_{k}= \Sigma_{i=0}^{N-1} f_{l}\omega_{N}^{-kl}');

参考：数式を使用したプロットへのラベルおよび注釈の付加　LaTeX を使用した数式を含むテキスト

2016年10月29日土曜日

Octaveレッスン(4) - 正弦波を生成する関数

* 正弦波を生成する関数

指定したサンプリングレート＆周波数の正弦波を生成し、ファイルに保存する

---- gen_sin.m の中身 ----

function[y, t] = gen_sin(Fs, hertz, seconds)
if(nargin == 0)
disp("usage: gen_sin(Fs, hertz(, seconds))");
return
elseif(nargin < 3)
seconds=1;
fprintf('default seconds = %d\n', seconds)
end

t=[0:seconds*Fs-1]*1/Fs; %time vector
y=sin(2*pi*hertz*t);
file=input('output file:')
if(sizeof(file)==0)
disp('not output')
else
audiowrite(file, y, Fs)
end
end

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

実行例1：サンプリングレート44.1k, 10Hzの正弦波を2秒生成し、ファイル保存し、描画

>> [y, t] = gen_sin(44100, 10, 2);
output file:'sin660Hz.wav'
file = sin660Hz.wav
>> plot(t,y)

実行例2：サンプリングレート44.1k, 2Hzの正弦波を2秒生成し、前のグラフに重ねて描画

>> [y2, t] = gen_sin(44100, 2, 2);
output file:
file = [](0x0)
not output file
>> hold on
>> plot(t,y2,'r-')
>> hold off

グラフ表示ウィンドウ上部の[Z+] でズーム、選択した領域をズームなど可能。

参考：波形生成: 時間ベクトルと正弦波　オーディオファイルの書き込み

2016年10月28日金曜日

Octaveレッスン(3) - wavファイルを入力しFFTをかける関数（サブ関数使用）

* wavファイルを入力しFFTをかける関数（サブ関数使用）

ファイル先頭の基本関数の下には、サブ関数を定義可能。
サブ関数は、ファイル内の基本関数＆他のサブ関数からのみ呼び出し可能。

※Octaveレッスン(1), (2) の記述では、wavファイルが1秒以上の場合に対応していなかった。
　今回のsample_func2関数を使用すると、正しい結果が得られる。

---- sample_func2.mの中身 ----

function[Y, y, Fs] = sample_func2(wavefile, targetHeltz)
disp('sample_func2() returns FFT spectrum & sample rate')
info = audioinfo(wavefile)
[y,Fs] = audioread(wavefile);

Y=fft(y);

L = info.TotalSamples %sample length
T = 1/Fs; %sampling period
t = (0:L-1)*T; %time vector

P2 = abs(Y/L); %both side spectrum
P1 = P2(1:L/2+1); %one side spectrum
P1(2:end-1) = 2*P1(2:end-1); %why from 2?
f=Fs*(0:L/2)/L;

plot_spectrum(f, P1)
plot_magnified(f, P1, Fs, L, targetHeltz-100, targetHeltz+100)
end

function[] = plot_spectrum(f, P1)
subplot(2,1,1)
title('spectrum')
xlabel('(Hz)')
plot(f,P1)
end

function[] = plot_magnified(f, P1, Fs, L, minf, maxf)
subplot(2,1,2)
title('magnified spectrum')
xlabel('(Hz)')
grid('on')
min=L/Fs*minf;
max=L/Fs*maxf;
plot(f(min:max),P1(min:max))
end

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

参考：スクリプトと関数

登録: 投稿 (Atom)