>> 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')
>> pkg install -forge ltfatドキュメントを読む限り、Mac版もWindows版も同じ手順でインストールできないとおかしいのだが、別のMacや旧いOSバージョン、旧いOctaveバージョンでも同様に失敗する。
>> ltfatstartと実行すれば、各種 ltfat の関数が利用可能になる。
>> demo_waveletsさらに処理速度を高速化したい場合には、
>> ltfatmexと実行すればよい・・・ ハズなのだが失敗する。
$ ./configure --enable-threads
$ make
$ make clean
$ ./configure --enable-float --enable-threads
$ make
>> ltfatmex
>> pkg load signalMacの場合は、まずinstallする必要がある：
>> pkg install -forge signal上記は、ネットワーク経由でダウンロード＆インストールするコマンドだが、おそらくcontrol package が未インストールのため失敗するだろう。
>> pkg install -forge control依存関係のあるcontrolをインストールできたので、signalをインストールしてロードする：
>> pkg install -forge signal
>> pkg load signalこれで、xcorr関数などが利用可能になる。
>> pkg listまた、別 package 無しで利用可能な関数の一覧は、[ドキュメント]タブ - 「Function Index」sectionで確認可能。
[pks, locs] = findpeaks(S,"DoubleSided", "MinPeakHeight", 0.3);と実行すると、pksに±0.3以上の極値（極大値＆極小値）のリストが返され、locsに極値のインデクスが返される。
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
F(ω)=∫f(t)e^(-jωt) dt …(1)上記2式は下式のように略記されることが多い：
f(t)=1/2π∫F(ω)e^(jωt)dt …(2)
F(ω)=F{f(t)}線形システムにおいて、インパルス応答h(t)のフーリエ変換をH(ω)とすると、下式が成り立つ：
f(t)=F^-1{F(ω)}
f(t)＊h(t)=1/2π∫F(ω)H(ω)e^(jωt)dt時間領域の畳み込み演算（＊）だが、周波数領域では単なる積の演算になっている。
h(t)=1/πt sin(ωt)実用の場では、離散フーリエ変換を使用するので、下図の右上のような離散的なsinc関数になる。
2πf(-ω)=∫F(t)e^(-jωt)dtつまり、フーリエ変換の対称性を示す下式が得られた：
= F{F(t)}
F{F(t)} = 2πf(-ω)フーリエ変換の対称性が役立つのは、f(t)のフーリエ変換F(ω)が判明しているが、F(ω)のωをtに入れ替えた関数F(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(ω) = 2πr(-ω)/2 = πr(-ω)と簡単に求まる。
pkg load signal次数10、低域カットオフ周波数 440 Hz 、高域カットオフ周波数 660 Hzのバンドパスフィルターを作成する例：
Fs = 44100; % sampling rate[b, a] に、伝達関数の係数が返される。
[b, a] = butter(10, [440 660]/(Fs/2));
[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
[y, t] = gen_sin(Fs, hertz, seconds);実行すると、wavファイル作成には失敗したが（audiowrite関数がない）、出力 y, t を得られた。
 |
| 『これなら分かる応用数学教室』 p.130 |
>> 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}');
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
>> [y, t] = gen_sin(44100, 10, 2);
output file:'sin660Hz.wav'
file = sin660Hz.wav
>> plot(t,y)
>> [y2, t] = gen_sin(44100, 2, 2);
output file:
file = [](0x0)
not output file
>> hold on
>> plot(t,y2,'r-')
>> hold off
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
function [出力変数] = 関数名(入力変数)ファイル名は関数名と合わせることが推奨されている。
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
>> file='sin660.wav'
>> Y=sample_func(file);
>> 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
>> L = info.SampleRate
>> T = 1/Fs
>> t = (0:L-1)*T
>> plot(t, y)
>> Y=fft(y)
>> 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
>> plot(f,P1)
>> plot(f(1:500),P1(1:500))
