信号処理のお仕事メモ: 10月 2016

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

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

参考：スクリプトと関数

Octaveレッスン(2)

* ファイル(*.m)に関数を定義して実行する

wavファイルを読み込んで、FFTをかけてグラフ表示する関数の例。

一行目の関数宣言には、複数の出力変数・入力変数を指定可能：

function [出力変数] = 関数名(入力変数)

ファイル名は関数名と合わせることが推奨されている。

---- sample_func.mの中身 ----

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

Y=fft(y);

L = info.SampleRate
T = 1/Fs;
t = (0:L-1)*T;

P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f=Fs*(0:L/2)/L;

subplot(2,1,1)
title('spectrum')
xlabel('(Hz)')
plot(f,P1)

subplot(2,1,2)
title('magnified 500Hz-700Hz')
xlabel('(Hz)')
grid on
plot(f(500:700),P1(500:700))
end

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

* .mファイルに定義した関数を実行する

sample_func.mで定義したsample_func関数を実行する例。
ローカルディレクトリにある.mファイルが自動的に読み込まれる。

>> file='sin660.wav'
>> Y=sample_func(file);

2016年10月27日木曜日

Octaveレッスン(1)

* Audacityで正弦波の音声ファイル(wav)を作成する

正弦波の波形を生成： [ジェネレーター] - [トーン…] で、Hz数、長さなどを設定
音声を確認する（繰り返し再生）： space + return押し → space押しで停止
WAVファイルとして保存： [ファイル] - [オーディオの書き出し] - wavファイル形式を指定

* Octaveでwavファイルを読み込む

>> pwd
ans = /Users/MyName
>> ls
sin440.wav sin8800.wav
>> info = audioinfo('sin440.wav')
info =
scalar structure containing the fields:
Filename = sin440.wav
CompressionMethod =
NumChannels = 1
SampleRate = 44100
TotalSamples = 44100
Duration = 1
BitsPerSample = 16
BitRate = -1
Title =
Artist =
Comment =
>> [y, Fs] = audioread('sin440.wav'); % y=samples, Fs=sample rate
>> sound(y,Fs) % play sound