息抜きで、以前の図形を、GIFアニメーションで連結。
諸事情により途中で終了。
「Matlab」タグアーカイブ
Matlab で マンデルブロ集合
息抜きで、Matlabでマンデルブロ集合の図形を描写するプログラムを作成してみました。
マンデルブロ集合とは
マンデルブロ集合は、Wikipediaから以下です。
次の漸化式で定義される複素数列
が
の極限で無限大に発散しないという条件を満たす複素数
全体が作る集合がマンデルブロ集合である。
（Wikipediaより引用）
![\[ \left\{ \begin{aligned} z_{n+1} &= {z_n}^2 + c \\ z_0 &= 0 \end{aligned} \right. \]](https://www.lbtora.net/wp-content/ql-cache/quicklatex.com-e7db5815dfee20b84c0352ae749c2865_l3.png "Rendered by QuickLaTeX.com")
ざっくりとしたプログラムの解説としては、複素平面座標にて、その複素数の発散の度合いで図形が描ける、というものです。
ソースコード
実際のソースコードは以下です。
複素平面にて描写範囲を指定し、各座標の複素数
の発散を調べています。
絶対値が2以上になった場合を発散しているとし、その発散回数を保存しています。
図形の描写では、各座標の発散までの計算回数をプロットしています。
（変数は単精度で計算しています）
% ----------------------------------------------------- %
% パラメータ設定
% ----------------------------------------------------- %
% 最大ループ回数
loop = 1e2;
% 各軸の描写点数
n = 1e3;
% 描写範囲
xmin = -1.5;
xmax = 0.5;
ymin = -1.0;
ymax = 1.0;
% ----------------------------------------------------- %
% ----------------------------------------------------- %
% マンデルブロ集合の計算
% ----------------------------------------------------- %
% 座標設定
x = single(linspace(xmin,xmax,n));
y = single(linspace(ymin,ymax,n));
C = ones(n,1)*x + 1i*y(:)*ones(1,n);
% マンデルブロ図形の計算
Z = zeros(n,'single');
W = zeros(n,'single');
% インデックスの初期値
index = true(size(C));
% ループ
for k = 1:loop
Z(index) = Z(index).^2 + C(index);
index = abs(Z) < 2;
W(~index & W==0) = k;
end
W(W==0) = loop+1;
% 図形のプロット
imagesc(x,y,W);
axis square;
title('Mandelbrot Set');
xlabel('Real');
ylabel('Imag');
% PNG 画像の保存
print('Mandelbrot.png', '-dpng');
% ----------------------------------------------------- %
このソースコードで描いてみた図形が以下。
描写位置を変更すると面白いです。
- Matlabのバージョンは 2014a
- マンデルブロ集合 (Wikipedia)
MatlabでLPF
備忘録です。
Matlabでローパスフィルタを簡易的に作成してみました。
- 処理手順(おそらく)
- LPFの変数指定
- FIRフィルタの係数計算
- フィルタ係数と入力信号から、通過信号算出
fdesign関数のパラメータをいじれば、いろいろなフィルタが簡易的にすぐに設計できそうです。
帯域を狭帯域にしすぎたり、阻止帯域での減衰が大きすぎると係数が増加し、処理時間も増加します。
function y = lpf_filtered_signal(x,Fs)
% LPF 処理信号算出関数
%
% 変数
% Fs : サンプリング周波数
% Fp : 通過帯域周波数
% Fst : 阻止帯域周波数
% Ap : 通過帯域のリップル [dB]
% Ast : 阻止帯域の減衰量 [dB]
% 変数の設定
Fp = 10*1e6;
Fst = 20*1e6;
Ap = 1;
Ast = 50;
h = fdesign.lowpass('Fp,Fst,Ap,Ast',Fp,Fst,Ap,Ast,Fs);
% フィルタの設計
H = design(h);
% フィルタ通過信号
y = filter(H.Numerator,1,x);
end
- 参考にしたサイト