高周波&数値計算関係記事リンク集

・Pythonの高周波用ライブラリ scikit-rfを使う

・Visual C#で数値計算ライブラリ math.NET numericsを使う(Visual Basic版もあり)

・高周波エンジニアのためのAI・機械学習入門

・JavaScriptの数値計算ライブラリmath.jsを使う

・カシオの高精度計算サイトに投稿した自作式一覧

» 続きを読む

2024年12月 9日 (月)

Pythonと高周波ライブラリscikit-rfを使って同軸コネクタ→同軸・基板接続不連続(LとC)→伝送線路(マイクロストリップライン)のような構成の測定治具のSパラメータ、特性インピーダンスを出力する関数を作る。後で機械学習に使うための準備。

LCフィルタ、パッチアンテナの次は伝送線路(マイクロストリップライン)と同軸コネクタ(SMAなど)を含む測定治具をやってみよう。

uSimmics(旧名QucsStudio)で以前やってみたこんな構成をPythonとscikit-rfで実施するイメージ。

Measurementjigsckitrf01

さっそく関数を紹介すると(もうコメントに説明があるので詳しくは略)、


import numpy as np
import matplotlib.pyplot as plt
import skrf as rf
from skrf.media import Coaxial, MLine
rf.stylely()

def mesurement_jig(fstart, fstop, n, coaxial_length, cap, ind, width, length, height, er):
    """
    測定治具を模したSパラメータ及び基板の特性インピーダンスを返す関数

    測定治具は
    同軸コネクタ→シャントC→シリーズL
    →基板に形成されたマイクロストリップライン
    →シリーズL→シャントC→同軸コネクタの構成
    導体損・tanδは今回は固定されている。

    Parameters
    ----------
    fstart : float
        最低周波数[GHz]
    fstop :  float
        最高周波数[GHz]
    n : int
        周波数分割数
    coaxial_length : float
        同軸コネクタ部分の長さ[mm]
    cap : float
        同軸コネクタ・基板接続部キャパシタンス[pF]
    ind : float
        同軸コネクタ・基板接続部インダクタンス[nH]
    width : float
        線路幅[mm]
    length : float
        線路長さ[mm]
    height : float
        基板厚み[mm]
    er : float
        基板比誘電率

    Returns
    -------
    Stl : scikit-rfのNetwork
        測定治具のSパラメータ
    z0 : float
        中央周波数の特性インピーダンス
    """
    #周波数範囲設定
    freq = rf.Frequency(fstart, fstop, n, "GHz")

    #同軸コネクタのパラメータ(SMA相当)
    coax = Coaxial(frequency=freq, Dint=1.3e-3, Dout = 4.59e-3, epsilon_r=2.29, tan_delta=4e-4, sigma=1/0.022e-6, z0_port=50)

    #マイクロストリップラインのパラメータ
    msl = MLine(frequency=freq, z0_port=50, w=width*1e-3, h=height*1e-3, t=35e-6, ep_r=er, tand=0.01, rho=1e-8, rough=0.127e-6)

    #同軸コネクタの長さ決定
    coax_line = coax.line(coaxial_length, unit="mm", name="coax_line")

    #マイクロストリップラインの長さ決定
    msl_line = msl.line(length,  unit="mm", name = "msl_line")

    #同軸コネクタと基板の接続部のLC
    C = msl.shunt_capacitor(cap * 1e-12)
    L = msl.inductor(ind * 1e-9)

    #Casccade接続する
    Stl = coax_line ** C ** L ** msl_line ** L ** C ** coax_line
    Stl.name = "Measurement jig"

    #中心周波数の基板の特性インピーダンスを求める。
    z0 = msl.z0[n // 2].real

    return Stl, z0

これで、例えば基板のインピーダンスは同じとしてL,Cを変えた場合のS11を見ると、

Stl1, z0_1 = mesurement_jig(0.1, 20, 200, 10, 0.05, 0.2, 0.49, 100, 0.254, 4)
Stl2, z0_2 = mesurement_jig(0.1, 20, 200, 10, 0.2, 0.2, 0.49, 100, 0.254, 4)
Stl3, z0_3 = mesurement_jig(0.1, 20, 200, 10, 0.05, 0.3, 0.49, 100, 0.254, 4)
Stl1.plot_s_db(m=0, n=0)
Stl2.plot_s_db(m=0, n=0)
Stl3.plot_s_db(m=0, n=0)
Measurementjigsckitrf02
とかなり違う。しかし基板部分の特性インピーダンスは
print(z0_1, z0_2, z0_3)
50.23630901433205 50.23630901433205 50.23630901433205
ともちろん全部一緒。
S11から特性インピーダンスを求めるにはTDR(Time Domain Reflectometry)を使うことが多いが、これからやろうとしていることはTDRを使わずにSパラメータだけから特性インピーダンスを機械学習で予測するということ。
では次回に。

高周波・RFニュース2024年12月9日 iFixitがDJI Neo分解、TechInsightsがApple Pencil Pro分解、QualcommのNeurIPS 2024でのAI技術発表、IntelのIEDM 2024での発表、 Nokiaの7GHz帯の6G、Analog DevicesのPhased Array Antennaのホワイトペーパー、ZDTが史上二番目の売上高

・iFixitがDJI Neo分解

DJI Neo Teardown: This Year’s Most Giftable Drone Has an Unfixable Camera

202412091

・TechInsightsがApple Pencil Pro分解

Apple Pencil Pro 2024 (A2538) Teardown

202412092

・QualcommのNeurIPS 2024でのAI技術発表

Qualcomm at NeurIPS 2024: Our groundbreaking innovations and cutting-edge advancements in AI

202412093

・IntelのIEDM 2024での発表

Intel Foundry Unveils Breakthroughs in Interconnect Scaling for Future Nodes

・Nokiaの7GHz帯の6G

Nokia kicks off 6G test in 7GHz at Dallas HQ

・Analog DevicesのPhased Array Antennaのホワイトペーパー

Highly Integrated Multibeam Beamformers Offer SWaP Benefits for Payload Phased Array Antennas

・ZDTが史上二番目の売上高

Zhen Ding Releases November 2024 Monthly Revenue Report

2024年12月 8日 (日)

松のやで極厚ロースかつ丼をいただく。確かにめちゃくちゃ厚くて、かつ熱くて食べ応え十分。

この写真ではわからないですが、めちゃくちゃ厚いです。それがめちゃくちゃ(ご飯も含め)熱々で出てきた。

これは食べ応えあった。三つ葉もいい感じ。

20241130-124702

2024年12月 7日 (土)

映画「ジャワーン(JAWAN)」を観てきた。もうなんでもありで面白かった!チャーリーズエンジェルからおやじの戦いからカーアクションから政治まで。タイトル回収と冒頭の子供の伏線回収にも痺れたし、敵役が悪すぎるのもカタルシスがある。もちろんダンスもあり。

ポストカードもらった。ちょっとネタバレしているような気もするが…

20241205-113654 20241205-162847

とにかく冒頭から「???」のシーンの連続で、30年後のハイジャックでさらに「???」となって途中である程度真相がわかるまではいったいこれはどういう話なんだと引き込まれる。

チーフと6人の女性(まさにチャーリーが出てくるチャーリーズエンジェルみたいな)の正体やその過去の壮絶さには涙…

とにかく悪役がめちゃくちゃ悪いです。よくこんな設定考えるなあというほど。

もちろんアクションもすごいですが、カーアクションまであるとは驚いた。これもすごかった。

でこれは公式サイトにも出てないのでネタバレ?なのかもですが主役のシャー・ルク・カーンさんの一人二役がびっくり。いちいち親父の登場がかっこいいが主役を食ってるというか本人だから…体もムキムキでかっこいい。

そしてタイトル回収シーンと、冒頭の伏線が回収されるうシーンは痺れた。

政治批判のようなこともあり、すごい勧善懲悪なのでインドで人気が出たのはよくわかる。

もちろんダンスもあり。面白いのでお勧め。

冒頭にdesclaimerが出るのと、エンドクレジットまったく名前読ませる気がないのには驚いたが…

2024年12月 6日 (金)

MATLAB Onlineで高周波基板設計用のRF PCB Toolboxを使ってみる。Coupled line バンドパスフィルタやratraceカプラが設計できる。モーメント法(MoM)や有限要素法(FEM)でちゃんと計算してくれているようだ。

RF Toolboxの次はRF PCB Toolbox。多層のP板の設計に使えるが、平面回路のデザインとソルバーも入っている。

まずはCoupled Line Filterをやってみる。チュートリアルから周波数だけ変えたもの。

まずはどんな構造か図示できる。

Matlabrfpcb003

そしてdesignで設計、Sパラメータに直す(結構時間かかる。ちゃんとシミュレーションしてくれているよう)。

Matlabrfpcb004

ラットレースカプラは基板をRogersに変えてみた。

Matlabrfpcb002

結果はこちら。ちゃんと3GHzで動くようになっている。

Matlabrfpcb001

ソルバーについてはこちら。

Overview of Solvers

しかしMATLABなんでもできるな。Antenna Toolboxもある。

 

 

高周波・RFニュース 2024年12月6日 NGMNが無線パフォーマンス評価フレームワーク発行、5GAAがC-V2Xのロードマップ発行、Marvellの3nm 1.6Tbps PAM4インターコネクト、Nokiaの2.4Tbps光伝送、Silicon Labsの低消費電力モジュール、Xiaomi 14T Pro分解動画

・NGMNが無線パフォーマンス評価フレームワーク発行

NGMN releases “Radio Performance Assessment Framework” to guide next-generation RAN development 

202412061

・5GAAがC-V2Xのロードマップ発行

5GAA publishes updated Roadmap for C-V2X

202412062

・Marvellの3nm 1.6Tbps PAM4インターコネクト

Marvell Unveils Industry’s First 3nm 1.6 Tbps PAM4 Interconnect Platform to Scale Accelerated Infrastructure 

・Nokiaの2.4Tbps光伝送

Nokia and Aramco successfully achieve first 2.4Tbps optical transmission

・Silicon Labsの低消費電力モジュール

Silicon Labs' Breakthrough Ultra-Low Power Wi-Fi 6 and Bluetooth LE 5.4 Modules Supercharge Device Deployment

202412063

・Xiaomi 14T Pro分解動画

 

2024年12月 5日 (木)

MATLAB Onlineで高周波用のRF Toolboxを使ってみる。Touchstoneファイルの読み込み、dB表示グラフ、スミスチャートなど簡単にできるし、フィルタ合成やIEEE P370 De-embedding(ZC-2xThru)も使える(MATLABで書かれたものがオリジナル)。

Interface誌でMATLAB Onlineのライセンスがついてきたのでいろいろ試している。

今回は高周波用途で使われるRF Toolbox。

読み込みは s = sparameters("ファイル名")で簡単にできるし、dB表示の図示はrfplot(s)でOK。ポートを指定する時は下図のようにする。

Matlabrf001

スミスチャートはsmithplot(s, i, j)とこれも簡単に描ける。

Matlabrf002

フィルタ合成もできる。

Matlabrf004

IEEE P370 De-embedding(ZC-2xThru)も使える。

Matlabrf003

もともとこれはMATLABで書かれていたので当然といえばそう。

Pythonのscikit-rfよりずっと簡単な気がする。ただTDRがどうも一回有理関数近似して、とか制限あるっぽい(自分でFFTで作ればいいだけですが)。RF pcb toolboxというのもあって、それを使えば複雑な高周波の平面回路も簡単にできそうなのでそれも試したい。

 

 

 

 

2024年12月 4日 (水)

高周波回路シミュレータQucsStudioがuSimmicsに名称変更し、バージョンも4.8.3から5.8にアップデートされた。Qucsと区別するためだそうだ。また、Pythonの高周波用ライブラリscikit-rfもv1.5.0にバージョンアップされていた

数日前に、QucsStudioの最新版が出た、と思ったらuSimmicsに名称変更されていた。

https://qucsstudio.de/download/

作者さんによると、

https://qucsstudio.de/forums/topic/version-5-8-released/」

待望の新バージョンがダウンロード セクションで入手可能です。Qucs プロジェクトと明確に区​​別するために、アプリケーション名が uSimmics に変更されたことに注意してください。これは、非常に多くの改良、新機能、バグ修正を含むメジャー リリースです。

ということでした。

画面の比較をすると、

QucsStudio

Qucsstudio

uSimmics

Usimmics

ということでアイコンがでかくなった以外はそこまで見た目は変わらない。ただ新機能が追加されているということなのでまた調べてみよう。

また、Pythonの高周波用ライブラリscikit-rfもv1.5.0にバージョンアップされていた。

https://github.com/scikit-rf/scikit-rf/releases/tag/v1.5.0

こちらもpipでアップデート済み。

 

MATLAB OnlineのSimulinkでローレンツ方程式をode8で計算してみる。Interface 2025年1月号でMATLAB Onlineの半年ライセンスがついてきたので。Simulinkを使うのは初めてだったが、わかりやすいSimulink入門コースを修了したのですぐできた。

さて昨日はInterfaceを買ってMATLAB Onlineのライセンスをゲットした話を書いた。

Interface2025年1月号はMATLABで1ニューロンから手作り 数学&図解でディープ・ラーニング。初歩からAlexNetの転移学習、CNNまで話題が豊富でなんとMatlab Onlineの半年ライセンスがついてくる。Simulinkや各種toolboxも使える。早速MATLAB入門オンラインコース修了した。

今日はその中でSimulinkを使ってLorenz方程式を解いてみよう。実はMATLABは昔使っていたが、Simulinkを使うのは今回が初めて。

オンライン講座のSimulink入門がよくできていたので簡単にできた。

ローレンツ方程式は

dx/dt = σ(y - x)

dy/dt = x(ρ - z) - y

dz/dt = xy - βz

の形をしていて、それをSimulinkに落とし込むとこうなった。

Simulinklorenz02

ソルバーは指定できるのでode8にした。結果はこちら。

Simulinklorenz01

箱を線で結ぶだけでできるのは(今更ながら)新鮮。

データインスペクターを使うと重ねて描ける。

Simulinklorenz03

ついでにロジスティック写像の分岐図も描いてみた。

 

2024年12月 3日 (火)

日経サイエンス2025年1月号の特集 和算再発見の佐藤賢一さんの記事「算聖 関孝和の実像」に出てきた矢高に対する円弧の2乗の近似式をカシオの高精度計算サイトkeisan.casio.jpの自作式として作った。ものすごい精度であることがよくわかる。

日経サイエンス2025年1月号の特集は和算再発見だった。

20241203-174049

アマゾンリンク:https://amzn.to/41cyaa9

その中の佐藤賢一さんの記事「算聖 関孝和の実像」に出てきた円弧の2乗の矢高に対する近似式、

関・建部の研幾算法での方法
F(x)=(-9/40+6131089x/120+752638x²/45+340127x³/30-5966x⁴/3+699904x⁵/45)/113^2
関の括要算法での方法
H(x)=(5107600x-23835413x²+43470240x³-37997429x⁴+15047062x⁵-1501025x⁶-281290x⁷)
/(1276900*(1-x)⁵)
が正確な値
cos⁻¹(1-2x)²

と比べてあまりにもすごいのでカシオの高精度計算サイトkeisan.casio.jpに自作式として作ってみた。

リンクはこちら。

 関孝和と建部賢弘の和算での矢高と弦長から弧長を求める精度

こんな感じの画面になります。

Wasan03

Wasan01

Wasan02

やっぱり関孝和はすごいわ。

Interface2025年1月号はMATLABで1ニューロンから手作り 数学&図解でディープ・ラーニング。初歩からAlexNetの転移学習、CNNまで話題が豊富で、なんとMatlab Onlineの半年ライセンスがついてくる。Simulinkや各種toolboxも使える。早速MATLAB入門オンラインコース修了した。

今月のInterfaceはMATLABで1ニューロンから手作り 数学&図解でディープ・ラーニング。

20241202-174726

アマゾンリンク:https://amzn.to/3OzDahr

Simulinkを使って本当に1ニューロンからモデル化してとても分かりやすい。

Matlabintroduction02

そこからAlexNetの転移学習や、CNNの実装、バックプロパゲーション、ADAM、また数学の基礎まで内容が豊富。しかしPythonで私もディープラーニングやっているけどMATLABのほうがはるかに簡単に記述できるな。

そして何と、2025年5月までのMATLAB Onlineのライセンスがついてくる!それ以外にも

・MATLAB Online ・Simulink
・Deep Learning Toolbox
・Statistics and Machine Learning Toolbox
・Image Processing Toolbox
・Computer Vision Toolbox
・Symbolic Math Toolbox

が使える。しかしMATLABは大昔使ったっきりで完全に忘れているのでオンラインコース受けてみた。

MATLAB入門のコース修了証。

Matlabintroduction01

Simulink入門の修了したのだが、うまくコース修了証が表示されない…

ディープラーニング以外にもいろいろ遊んでみよう。

知らなかったのだが、MATLAB Onlineはライセンスなしでも無償で月20時間使えるそうだし。

また付録はディープ・ラーニングの始まりと現代社会での活用、だったがこれも読みごたえがある。エミー・ネーターまで記載されていたのは驚いた。

2024年12月 2日 (月)

高周波エンジニアのためのAI・機械学習入門(18)パッチアンテナ設計モジュールをPythonで作ったので、電気的に測定できる入力インピーダンス、共振周波数を与えて指向性(Directivity)をKeras3.0を用いたディープラーニング(DNN)による回帰で予測する。

さて前回はパッチサイズ(W、L)、基板厚みh、共振周波数を与えて基板誘電率が予測できるかやってみた。

今回は電気的に測定できる入力インピーダンスRinと共振周波数f0のみを与えて指向性が予測できるかやってみる。

学習データはいつものようにモンテカルロシミュレーションで作る。


import numpy as np
import patchantenna as pa

#モンテカルロシミュレーションでは基板比誘電率、厚み、中心周波数を乱数で振るが、それの最大最小
ermin = 2.0
ermax = 10.0
hmin = 0.1
hmax = 1.0
f0min = 1
f0max = 10

#モンテカルロシミュレーションでデータ作成して保存する。
#周波数と入力インピーダンスを与えて指向性を予測する。
N = 10000
np.random.seed(1)
data = np.empty((0, 2))
label = np.empty((0, 1))
for i in range(N):
    er = ermin + (ermax - ermin) * np.random.rand()
    h = hmin + (hmax - hmin) * np.random.rand()
    f0 = f0min + (f0max - f0min) * np.random.rand()    
    ant = pa.PatchAntenna(er, h, f0)
    Rin = ant.input_impedance()
    direc = ant.directivity()
    data = np.vstack((data, np.array([Rin, f0])))
    label = np.vstack((label, direc))

#保存
np.savez_compressed("directivity.npz", data=data, label=label)

 

そしてKeras3.0を用いたディープラーニング(DNN)で予測する。
まずデータを読み込み、訓練データとテストデータに分ける。


import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
import os
os.environ["KERAS_BACKEND"] = "tensorflow"
import keras

data_label = np.load("directivity.npz")
data = data_label["data"]
label = data_label["label"]
x_train, x_test, y_train, y_test = train_test_split(data, label, test_size=0.3, random_state=0)

次にモデルを作る。セル数64の3層のモデル。


#正規化
normalizer = keras.layers.Normalization()
normalizer.adapt(x_train)

# Functional APIでDense層を3層にしたDNNを設定
hidden_dim = 64
inputs = keras.Input(shape=(2,))
x = normalizer(inputs)
x = keras.layers.Dense(hidden_dim, activation="relu")(x)
x = keras.layers.Dense(hidden_dim, activation="relu")(x)
x = keras.layers.Dense(hidden_dim, activation="relu")(x)
outputs = keras.layers.Dense(1)(x)

# モデルの設定
model = keras.Model(inputs=inputs, outputs=outputs)
model.compile(loss = 'mean_squared_error' ,optimizer=keras.optimizers.Adam())

そして学習させる。


batch_size = 32
epochs = 100

keras.utils.set_random_seed(1)
history = model.fit(
    x_train,
    y_train,
    batch_size=batch_size,
    epochs=epochs,
    validation_split=0.15,
)
y_pred = model.predict(x_test)
metric = keras.metrics.R2Score()
metric.update_state(y_test, y_pred)
result = metric.result()
print(result)
error = np.abs((y_test - y_pred)/y_test*100)
print(error.mean(axis=0))

R2score99.7%となかなかいい結果になった。

では図示して確認してみる。


legend = ["Directivity"]
fig, ax = plt.subplots(1, 2, figsize=(12,6))

maxvalue = y_pred.max()
ax[0].scatter(y_pred, y_test, c="r", s=5)
ax[0].plot([0,maxvalue], [0,maxvalue], "--", c="black")
ax[0].set_xlabel("推定した値", fontname='MS Gothic')
ax[0].set_ylabel("実際の値", fontname='MS Gothic')
ax[0].set_xlim(5, maxvalue)
ax[0].set_ylim(5, maxvalue)
ax[0].grid()
ax[0].legend([legend[0] + f" 平均誤差{error.mean():.2f}%"], prop={"family":"MS Gothic"})
ax[1].hist(error, bins = 100)
ax[1].set_xlabel("誤差[%]", fontname='MS Gothic')
ax[1].set_ylabel("頻度", fontname='MS Gothic')
ax[1].grid()
fig.tight_layout()
plt.show()

 

Dnndirectivity

これだけ合っていれば使えそう。次の題材は伝送線路かな。

 

2024年12月 1日 (日)

神座で麻婆麵をいただく。しっかり花椒が効いていて、辛さもありなかなか美味しかった。もちろん白菜もたっぷり。

前々から気になっていたがようやく神座で麻婆麵を食べた。

20241129-112442

想像していたよりずっと本格的な麻婆豆腐が白菜たっぷりのおいしいラーメンの上にかかっている感じ。

ちゃんと花椒も感じられるし、辛さもある。肉も多い。これは当たりだった。

 

2024年11月30日 (土)

「アングリースクワッド 公務員と7人の詐欺師」を観てきた。スカッとするお話で面白かった!内野聖陽さんが全編顔芸ですごかったし、岡田将生くんの悪そうな顔も小澤征悦さんの悪役もよかった。で途中で数えてあれ?と思ったり、最後にえ?と思ったり。さすがはカメ止めの上田監督。

有名俳優さんが惜しげもなく使われていて、上田監督もすごくなった(上から目線で済みません…)とまず最初に思ったり。

あとパルムが食べたくなる。

20241125-123731

とにかく内野聖陽さん、ずっと顔芸しててめちゃくちゃ面白かった。岡田将生くんも、完全な二枚目よりはこういうちょっと悪いところがある役は合ってると思った。

でとにかく小澤征悦さんの悪役と吹越満さんの小悪党が悪くて悪くて、これは最後にどうやってだます?と思ったら大ピンチでドキドキ。

そこはさすがの上田監督であっと驚く展開でスカッとします。

最後にえ?というところがあるんですが、あれはいるのかな…本当の最後だけでよかった気も…でもあれで逆にスカッとする人もいるだろうからありなのかな。

まあとにかく面白いのでお勧めです。

2024年11月29日 (金)

線表現の可能性@国立国際美術館を観てきた。子供が描いたような線、緻密な線(ジョナサン・ボロフスキーと池田龍雄がよかった)、カラフルな線、立体的な線(あれどうやってバランス?)など様々で面白かった。彼女の肖像も一緒に観られたがこれも興味深かった。

久しぶりに国立国際美術館へ。線表現の可能性を観に来た。

20241120-125540

 

20241120-125542

20241120-125911

写真撮影はOKでした。子供の落書きのような線から始まり、緻密だったり(ジョナサン・ボロフスキーのノートの切れ端に書いた線や、池田龍雄の化物の系譜が気持ち悪くてよかった)、カラフルなもの(ポスターにもなってるベルナール・フリズのガルプなど)、立体的なもの(植松奎二の置ー浮くかたちがどうやってバランスを取っていたのか全然分からん…)など様々で面白かった。
後半は2020年代に亡くなった方々の作品で、イリヤ・カバコフの天使と出会う方法がよかった。

その次はコレクション1 彼女の肖像へ。

20241120-134710

岡本信治郎の版画集『ベティ・ブープの国』や小沢剛のベジタブル・ウェポンが面白かった。そして谷原菜摘子のSADOにはおおっと思った。人魚か。

ところで来ていたのはほぼ外国の方々だった。さすが国際美術館。

«高周波・RFニュース 2024年11月29日 NuvotronicsのPolystrataがQorvoのGaNパワーアンプに採用、WürthがBalunのラインアップ拡張、MinewSemiのモジュールにNordicのnRF54Lが採用、AIチップにPCIe 7.0が必要な理由、SIJの特集はThe Road from 1 Gbps-NRZ to 224 Gbps-PAM4

最近の記事

最近のコメント

2024年12月
1 2 3 4 5 6 7
8 9 10 11 12 13 14
15 16 17 18 19 20 21
22 23 24 25 26 27 28
29 30 31        
フォト
無料ブログはココログ