2025-07-06

Plotlyでエクセル上のデータをプロットするスクリプトを作成した。

はじめに

Excel は表計算やグラフ作成、数式によるデータ演算、条件付き書式を使ったヒートマップ生成など、非常に便利なツールである。
しかし個人的には、グラフ描画機能に以下のような非効率な側面を感じている。

Excel のグラフ描画機能における課題

1. プロットデータの登録方法が非効率的

GUI 上でセルを１つずつ選択しながらプロット範囲を登録する必要がある
「グラフを右クリック → データの選択 → 挿入 → 範囲指定」と、多くのマウス操作が発生
データが連続していればまだしも、飛び飛びに入力されたセルを登録するのは非常に手間がかかる

2. デフォルトテンプレートの設定が変更できない

定義済みのテンプレートを呼び出すことは可能だが、既定のグラフスタイル（軸ラベルの位置や枠線の太さなど）を変更できない
理工系には独自の「ローカルルール」があってもおかしくないが、Excel では毎回手動で調整せざるを得ない

3. 描画範囲の調整に手間がかかる

MATLAB や Python のインタラクティブプロット、LabVIEW のように、マウスドラッグでズームできない
「グラフの書式設定」ダイアログから最小値・最大値を手入力する必要があり、ちょっと波形を拡大したいだけでも二手間以上かかってしまう

上記のような不満は、理工系分野で Excel を多用してきた方なら共感いただけるのではないだろうか。

Python＋Plotly で散布図を自動生成するスクリプト

今回は、インタラクティブかつ美しい体裁のグラフが簡単に作れる“Plotly”ライブラリを用いて、
Excel に入力されたデータを半自動でプロットする Python スクリプトを紹介する。jupyter notebookで利用するとより便利でおススメである。

セル範囲はコード中に文字列で指定
完全自動化ではないが、データ範囲が規則的であれば「繰り返し実行」→「再現性担保」→「効率化」につながる

本スクリプトで上記の課題がすべて解決できるわけではないが、
データ処理やグラフ作成の作業効率を向上させる一助となれば幸いである。

使い方例1

使い方はいたってシンプルである。

① 関数generate_plot_infoの引数にシート名、凡例セル、x軸の列名、y軸の列名、データ行の範囲を渡して、それをplot_info_list配列に格納する。以下の例ではプロットするデータが一つのため、plot_info_list配列には一つしかデータが登録されていないが、複数データをプロットする場合は、プロットするデータ数に合わせてgenerate_plot_info関数の返却値をplot_info_list配列に追加する。たとえば以下のようなイメージ。

plot_info_list = [
    generate_plot_info(sheet_name=data_com_info["sheet_name"], legend_cell=data_com_info["legend_cell"], x_col=data_com_info["x_col"], y_col=data_com_info["y_col"], data_range=data_com_info["data_range"]),
    generate_plot_info(sheet_name=data_com_info["sheet_name"], legend_cell=data_com_info["legend_cell"], x_col=data_com_info["x_col"], y_col=data_com_info["y_col"], data_range=data_com_info["data_range"]),
    generate_plot_info(sheet_name=data_com_info["sheet_name"], legend_cell=data_com_info["legend_cell"], x_col=data_com_info["x_col"], y_col=data_com_info["y_col"], data_range=data_com_info["data_range"])
]

② 関数plot_excel_dataの引数にファイルパス、前の手順で作成したデータ情報配列plot_info_list配列、軸範囲設定モード、グラフタイトル、x軸ラベル、y軸ラベル、を渡し、関数を実行する。
軸範囲設定モードについて補足する。モードに応じて、x,y軸のデフォルト範囲設定のどこを自動設定するかを設定できる。
設定値と設定内容の対応は以下の通り。
1: x,y ともに [-cv_max, cv_max]
2: x,y ともに自動
3: x 自動, y は y_custom
4: y 自動, x は x_custom
5: x,y ともに x_custom, y_custom

# Example usage
# Excelファイルのパスとプロット情報を指定 
file_path = r"C:\Users\User\Documents\test_waveform.xlsx"
title_info = {
    "graph_title": "Sample Plot",
    "graph_xaxis": "X Axis",
    "graph_yaxis": "Y Axis"
}
data_com_info = {
    "sheet_name": "DRV0022_uvw_wvu0",  # シート名
    "legend_cell": "B1",  # 凡例セル
    "x_col": "F",  # x軸の列
    "y_col": "G",  # y軸の列
    "data_range": [2, 1001]  # データ範囲 [開始行, 終了行]
}
# プロット情報を生成
plot_info_list = [
    generate_plot_info(sheet_name=data_com_info["sheet_name"], legend_cell=data_com_info["legend_cell"], x_col=data_com_info["x_col"], y_col=data_com_info["y_col"], data_range=data_com_info["data_range"])
]
display(plot_info_list)
# プロット実行
plot_excel_data(file_path, plot_info_list, axis_mode=1, graph_title=title_info["graph_title"], graph_xaxis=title_info["graph_xaxis"], graph_yaxis=title_info["graph_yaxis"])

使い方例②：追加データを重ねてプロット

重ねてプロットの一例（青色データがExcelから抽出したデータ、赤色データはpythonコード上から作成したデータ）

# Example usage
# Excelファイルのパスとプロット情報を指定 
file_path = r"C:\Users\User\Documents\test_waveform.xlsx"
title_info = {
    "graph_title": "Sample Plot",
    "graph_xaxis": "X Axis",
    "graph_yaxis": "Y Axis"
}
data_com_info = {
    "sheet_name": "DRV0022_uvw_wvu0",  # シート名
    "legend_cell": "B1",  # 凡例セル
    "x_col": "F",  # x軸の列
    "y_col": "G",  # y軸の列
    "data_range": [2, 1001]  # データ範囲 [開始行, 終了行]
}
# プロット情報を生成
plot_info_list = [
    generate_plot_info(sheet_name=data_com_info["sheet_name"], legend_cell=data_com_info["legend_cell"], x_col=data_com_info["x_col"], y_col=data_com_info["y_col"], data_range=data_com_info["data_range"])
]
display(plot_info_list)
# プロット実行
fig = plot_excel_data(file_path, plot_info_list, axis_mode=1, graph_title=title_info["graph_title"], graph_xaxis=title_info["graph_xaxis"], graph_yaxis=title_info["graph_yaxis"], show_plot=False)
# 追加データのプロット
x_additional = np.linspace(-15, 15, 100)
y_additional = x_additional  # 例: y = x のデータ
fig.add_trace(
    go.Scatter(x=x_additional, y=y_additional, mode='lines', name='Additional Data')
)
fig.show()  # 最後に全てのデータを表示

関数スクリプト

import openpyxl
import plotly.graph_objects as go
import numpy as np

def plot_excel_data(
    file_path: str,
    plot_info_list: list,
    axis_mode: int = 1,
    graph_title: str = "Excel Data Plot",
    graph_xaxis: str = "X Axis",
    graph_yaxis: str = "Y Axis",
    x_custom: list[float] | None = None,
    y_custom: list[float] | None = None,
    c: float = 1.2,
    canvas_size: int = 600,
    show_plot: bool = True
):
    """
    Excel のデータを読み込んで Plotly でプロット。軸レンジは5モードから選択可。

    Args:
        file_path: Excel ファイルのパス
        plot_info_list: プロット情報のリスト。各要素は
            [sheet_name, legend_cell, x_range, y_range]
        axis_mode: 軸レンジモード (1〜5)
          1: x,y ともに [-c*v_max, c*v_max]
          2: x,y ともに自動
          3: x 自動, y は y_custom
          4: y 自動, x は x_custom
          5: x,y ともに x_custom, y_custom
        x_custom: モード4,5 で使う x 軸の [min, max]
        y_custom: モード3,5 で使う y 軸の [min, max]
        c: モード1 の倍率 (デフォルト 1.2)
        canvas_size: 図全体のピクセル幅・高さ (正方形)
    """
    wb = openpyxl.load_workbook(file_path, data_only=True)
    fig = go.Figure()
    all_vals: list[float] = []

    # 各プロット情報を処理
    for sheet, legend_cell, xrng, yrng in plot_info_list:
        if sheet not in wb.sheetnames:
            raise ValueError(f"シート '{sheet}' が見つかりません")
        ws = wb[sheet]
        legend = ws[legend_cell].value

        x_cells = ws[xrng]; y_cells = ws[yrng]
        x_data = [c.value for row in x_cells for c in row]
        y_data = [c.value for row in y_cells for c in row]
        all_vals += [abs(v) for v in x_data + y_data if isinstance(v, (int, float))]

        fig.add_trace(
            go.Scatter(x=x_data, y=y_data, mode='lines+markers', name=str(legend))
        )

    # 軸レンジ設定
    xaxis: dict = {}
    yaxis: dict = {}
    if axis_mode == 1:
        if not all_vals:
            raise ValueError("プロットできる数値データがありません")
        v_max = max(all_vals); lim = c * v_max
        xaxis["range"] = [-lim, lim]; yaxis["range"] = [-lim, lim]
    elif axis_mode == 3:
        if y_custom is None or len(y_custom) != 2:
            raise ValueError("mode 3 では y_custom に [min,max] を指定してください")
        yaxis["range"] = y_custom
    elif axis_mode == 4:
        if x_custom is None or len(x_custom) != 2:
            raise ValueError("mode 4 では x_custom に [min,max] を指定してください")
        xaxis["range"] = x_custom
    elif axis_mode == 5:
        if (x_custom is None or len(x_custom) != 2
            or y_custom is None or len(y_custom) != 2):
            raise ValueError("mode 5 では x_custom, y_custom ともに [min,max] を指定してください")
        xaxis["range"] = x_custom; yaxis["range"] = y_custom

    fig.update_layout(
        xaxis=xaxis,
        yaxis=yaxis,
        width=canvas_size,
        height=canvas_size,
        title=graph_title,
        xaxis_title=graph_xaxis,
        yaxis_title=graph_yaxis
    )

    if show_plot:
        fig.show()
    else:
        return fig

def generate_plot_info(
    sheet_name: str,
    legend_cell: str,
    x_col: str,
    y_col: str,
    data_range: list[int]
) -> list[str]:
    """
    プロット情報を生成するヘルパー関数。

    Returns:
        [sheet_name, legend_cell,
         f"{x_col}{開始行}:{x_col}{終了行}",
         f"{y_col}{開始行}:{y_col}{終了行}"]
    """
    start_row, end_row = data_range
    return [
        sheet_name,
        legend_cell,
        f"{x_col}{start_row}:{x_col}{end_row}",
        f"{y_col}{start_row}:{y_col}{end_row}"
    ]

以上、Python＋Plotly で Excel データを効率的に可視化するスクリプトの紹介である。ぜひお試しのうえ、業務効率化にお役立てください

2025-06-28

Pythonでタスク分類ツールを作った

タスク分類ツール

このツールは、タスクをジャンルとメトリクスで分類し、視覚的に評価するためのアプリケーションである。タスクを登録し、一覧表示や結果のグラフ化をすることができる。

github.com

機能

タスクの登録
タスクの一覧表示
タスクの編集
タスクのメトリクスによる評価
結果のグラフ表示
メトリクス
緊急度
労力
影響度
必要なエネルギー
モチベーション
自己犠牲感
ジャンル
仕事
家事
リフレッシュ
自己研鑽
交際関係
使用方法
タスクを登録するには、タスク名、ジャンル、メトリクスを入力し、「登録」ボタンをクリックする。
登録したタスクは一覧に表示される。タスクを選択して「編集」ボタンをクリックすると、内容を変更できる。
タスクのメトリクスによる評価を行い、結果をグラフで表示する。
開発環境
Python 3.x
PyQt5
pandas
plotly
インストール

必要なライブラリをインストールするには、以下のコマンドを実行する。

pip install PyQt5 pandas plotly

実行方法

python task_classifier_v2.py

2025-05-27

ノッチフィルタ（バンドストップフィルタ）の勉強備忘録

E. Mattingly氏(Harvard 大学)が公開しているOpen Source MPI装置(OPS-MPI)では、至る所でフィルタ回路が登場する。中でも今回は"ノッチフィルタ（バンドストップフィルタ）"について取り上げる。

下記はOPS-MPIのDrive回路図であり、破線で囲まれた右端のブロックは励磁コイルの等価回路を表す。中段の緑色破線で囲まれたブロックがノッチフィルタであり、励磁信号に粒子起因の高調波成分が含まれないように信号純度（励磁周波数f0とノイズの比率）を高めることを目的として配置されている。

OPS-MPIの励磁信号入力ブロックに用いられているフィルタ回路（文献[1]より引用）

バンドストップフィルタ自体は非常に単純で、LPFとHPFを並列に接続させることで、LPFとHPFの減衰領域の帯域のみ通さないという原理。低周波信号はLPF側をパスし、高周波信号はHPF側をパスし、狙いとする周波数帯域の信号は両方のパスからブロックされるというものである。

下記ページにはR, Cを用いた単純な構成例のノッチフィルタが紹介されている。
どうやら、ノッチしたい狙いの周波数を f_N とすると、

によりC_NとR_Nを定め、
・LPF側：R=R_N（そのまま）、C=2*C_N（2倍）
・HPF側：R=R_N/2（1/2倍、　C=C_N（そのまま）
とすることで構成出来るとのこと。

detail-infomation.com

論文中ではインダクタとキャパシタ用いたフィルタを構成することで、2f0, 3f0を効果的にフィルタできるようなノッチフィルタを構成しているとのことだが、LPFとHPFの並列組み合わせという基本的なところは同一のはずである。

今回は、勉強を兼ねてLTSpiceにてノッチフィルタのシミュレーションを実行してみた。
検証した回路は以下の回路図にて示した通り。最下段に示す回路がノッチフィルタであり、最上段と中段にはそれぞれノッチフィルタで用いているLPF, HPFをそれぞれ単独で抽出した回路を示している。ノッチフィルタの周波数特性が、LPFとHPFの減衰領域のかぶる部分に対応して落ち込む特性となっているのかどうかを確かめることを目的としている。

素子値は R = 1kΩ、C = 0.1µF としたため、ノッチ周波数 f_N = 1592 Hz となる。
下記表にはパラメータ計算結果を示すが、LPFのfc（カットオフ周波数）= 796 Hz、HPFのfc = 3183 Hz となっており、ちょうどノッチ周波数 f_N がその間の値となることが分かる。

ノッチフィルタ（最下段）、ノッチフィルタを構成するLPF（最上段）、HPF（中段）

以下はシミュレーション結果の周波数特性を示す。
緑（LPF）と青（HPF）の減衰領域がクロスする帯域において、赤（ノッチフィルタ）の振幅特性はちょうど落ち込むような特性となっており、特定帯域のみブロックできるフィルタを構成できていることを確認できた。

参考文献

[1] E. Mattingly, “Design, construction, and validation of magnetic particle imaging systems for rodent, primate, and human functional neuroimaging,” Ph.D. thesis, Massachusetts Institute of Technology, 2024

2025-05-16

【技術メモ】PythonでMPI Physics Demo を再現してみた-Part 2

makutsueeken5.hatenablog.com

続編
正弦波入力だけでなく、
・矩形波
・三角波
etc. といろいろ入力可能。

# File: applied_field.py
import numpy as np

mu_0 = 4 * np.pi * 1e-7  # permeability of free space

def generate_applied_field(omega_drive, Hamp_drive, Hdc_drive, T_cycles, Npts_cycles):
    """
    Generate time vector, applied AC field, and sampling frequency.
    Returns:
        t: time vector [s]
        Hac: applied field array [A/m]
        Fs: sampling frequency [Hz]
    """
    period = 2 * np.pi / omega_drive
    total_time = period * T_cycles
    t = np.linspace(0, total_time, Npts_cycles)
    Hac = Hamp_drive * np.sin(omega_drive * t) + Hdc_drive
    Fs = Npts_cycles / total_time
    return t, Hac, Fs

# File: square_wave.py
import numpy as np

def generate_square_wave(frequency, amplitude=1, offset=0, T_cycles=1, Npts_cycles=1000, duty=0.5, harmonics=10):
    """
    Generate a continuous-approximation square wave via Fourier series.

    frequency: fundamental frequency [Hz]
    amplitude: peak amplitude
    offset: DC offset
    T_cycles: number of cycles
    Npts_cycles: number of sample points
    duty: duty cycle (0-1)
    harmonics: number of odd harmonics to sum

    Returns:
        t: time vector [s]
        waveform: approximate square wave [same units as amplitude]
        Fs: sampling frequency [Hz]
    """
    period = 1.0 / frequency
    total_time = period * T_cycles
    t = np.linspace(0, total_time, Npts_cycles)
    # Fourier series sum of odd harmonics
    waveform = np.zeros_like(t)
    for n in range(1, 2*harmonics, 2):
        # sigma approximates the square wave
        p = 1
        sigma_n = (np.sin(np.pi * n /(2*harmonics) ) / (np.pi* n/(2*harmonics)))**p
        waveform += sigma_n * (1.0 / n) * np.sin(2 * np.pi * n * frequency * t)
    waveform = amplitude * (4.0 / np.pi) * waveform + offset
    Fs = Npts_cycles / total_time
    return t, waveform, Fs

# File: triangle_wave.py
import numpy as np
from scipy.signal import sawtooth

def generate_triangle_wave(frequency, amplitude=1, offset=0, T_cycles=1, Npts_cycles=1000):
    """
    Generate a triangular waveform using scipy.signal.sawtooth with width=0.5.
    """
    period = 1 / frequency
    total_time = period * T_cycles
    t = np.linspace(0, total_time, Npts_cycles)
    waveform = amplitude * sawtooth(2 * np.pi * frequency * t, width=0.5) + offset
    Fs = Npts_cycles / total_time
    return t, waveform, Fs

# File: chirp_signal.py
import numpy as np
from scipy.signal import chirp

def generate_chirp_signal(f0, f1, duration=1, Npts_cycles=1000, amplitude=1, offset=0, method='linear'):
    """
    Generate a frequency-swept chirp signal.
    f0: start frequency [Hz]
    f1: end frequency [Hz]
    duration: total time [s]
    Npts_cycles: number of sample points
    method: interpolation method ('linear', 'quadratic', etc.)
    """
    t = np.linspace(0, duration, Npts_cycles)
    waveform = amplitude * chirp(t, f0=f0, f1=f1, t1=duration, method=method) + offset
    Fs = Npts_cycles / duration
    return t, waveform, Fs

# File: morlet_wavelet.py
import numpy as np

def generate_morlet_wavelet(center_freq, sigma, duration=1, Npts=1000):
    """
    Generate a Morlet (Gabor) wavelet.
    center_freq: center frequency [Hz]
    sigma: Gaussian envelope standard deviation [s]
    duration: total duration [s]
    Npts: number of sample points
    """
    t = np.linspace(-duration/2, duration/2, Npts)
    wavelet = np.cos(2 * np.pi * center_freq * t) * np.exp(-t**2 / (2 * sigma**2))
    Fs = Npts / duration
    return t, wavelet, Fs

# File: mexican_hat_wavelet.py
import numpy as np

def generate_mexican_hat_wavelet(sigma, duration=1, Npts=1000):
    """
    Generate a Mexican Hat (Ricker) wavelet.
    sigma: scale parameter [s]
    """
    t = np.linspace(-duration/2, duration/2, Npts)
    factor = 2 / (np.sqrt(3 * sigma) * (np.pi**0.25))
    wavelet = factor * (1 - (t / sigma)**2) * np.exp(-t**2 / (2 * sigma**2))
    Fs = Npts / duration
    return t, wavelet, Fs

# File: haar_wavelet.py
import numpy as np

def generate_haar_wavelet(duration=1, Npts=1000):
    """
    Generate a simple Haar wavelet (step function).
    duration: total duration [s]
    """
    t = np.linspace(0, duration, Npts)
    wavelet = np.where(t < duration/2, 1, -1)
    Fs = Npts / duration
    return t, wavelet, Fs



# File: response.py
import numpy as np

# from applied_field import mu_0

def langevin(x, c=1, Ms=5e5, D=20e-9, kb=1.380649e-23, T=300):
    """
    Compute Langevin response for input field array x.
    c : concentration [mol/L]
    """
    Na = 6.022e23  # Avogadro's number
    V = np.pi * D**3 / 6
    m = Ms * V
    beta = mu_0 * m / (kb * T)
    coeff = (c * Na * 1000.0) * m
    bH = beta * x
    y = np.zeros_like(bH)
    for i, val in enumerate(bH):
        if val == 0:
            y[i] = 0
        else:
            y[i] = coeff * (1/np.tanh(val) - 1/val)
    return y


def generate_response(H_field, **kwargs):
    """
    Wrapper to generate magnetization response via Langevin function.
    """
    return langevin(H_field, **kwargs)


# File: fft_processing.py
import numpy as np

def compute_fft(signal, Fs):
    """
    Compute FFT of a time-domain signal.
    Returns:
        freq: frequency array [Hz]
        amplitude: amplitude spectrum
        phase: phase spectrum
    """
    N = len(signal)
    fft_vals = np.fft.fft(signal) * 2 / N
    freq = np.fft.fftfreq(N, d=1/Fs)
    amp = np.abs(fft_vals) * 2 / N
    amp[0] /= 2  # correct DC component
    phase = np.angle(fft_vals)
    return freq, amp, phase


# File: plotting.py
import numpy as np
import plotly.graph_objs as go
from plotly.subplots import make_subplots
# from applied_field import mu_0
# File: app.py
from dash import Dash, html, dcc
import numpy as np


def create_figure(t, Hac, H_resp, H_L, M_L, f_H, A_H, f_M, A_M):
    """
    Create a 2x2 subplot figure showing
    - M(H)
    - H_response vs t
    - H_applied vs t
    - Normalized FFTs
    """
    fig = make_subplots(
        rows=2, cols=2,
        subplot_titles=('M(H)', 'H_response(t)', 'H_applied(t)', 'FFT Amplitude vs Frequency'),
        vertical_spacing=0.1
    )
    # M(H)
    fig.add_trace(
        go.Scatter(x=H_L, y=M_L*mu_0, mode='lines', name='M(H)', line=dict(dash='dash')),
        row=1, col=1
    )
    # Sweep response
    fig.add_trace(
        go.Scatter(x=Hac, y=H_resp*mu_0, mode='lines', name='Response'),
        row=1, col=1
    )
    # H_applied(t)
    fig.add_trace(
        go.Scatter(y=t, x=Hac, mode='lines', name='H(t)'),
        row=2, col=1
    )
    # M(t)
    fig.add_trace(
        go.Scatter(x=t, y=H_resp*mu_0, mode='lines', name='M(t)'),
        row=1, col=2
    )
    # FFTs
    half = len(f_H)//2
    fig.add_trace(
        go.Scatter(x=f_M[:half], y=A_M[:half]/np.max(A_M[:half]), mode='lines', name='FFT(M)'),
        row=2, col=2
    )
    fig.add_trace(
        go.Scatter(x=f_H[:half], y=A_H[:half]/np.max(A_H[:half]), mode='lines', name='FFT(H)'),
        row=2, col=2
    )
    # Axis labels and layout
    fig.update_xaxes(title_text='H [T]', row=1, col=1)
    fig.update_yaxes(title_text='M [A/m]', row=1, col=1)
    fig.update_yaxes(title_text='t [s]', row=2, col=1)
    fig.update_xaxes(title_text='H [T]', row=2, col=1)
    fig.update_xaxes(title_text='t [s]', row=1, col=2)
    fig.update_yaxes(title_text='M [A/m]', row=1, col=2)
    fig.update_xaxes(title_text='Frequency [Hz]', row=2, col=2, type='log')
    fig.update_yaxes(title_text='Normalized Amp', row=2, col=2, type='log')
    fig.update_layout(height=800, width=1200, title_text='Langevin Function Response')
    # range for x-axis in subplot (2,1), (1,1), (1,2)
    max_scale = 1.2
    max_Hac = np.max(Hac)*max_scale
    max_Hresponse = np.max(H_resp*mu_0)*max_scale
    max_Mresponse = np.max(M_L*mu_0)*max_scale
    # xrange for H applied (2,1) is equal to xrange for M-H (1,1)
    fig.update_xaxes(range=[-max_Hac, max_Hac], row=1, col=1)
    fig.update_xaxes(range=[-max_Hac, max_Hac], row=2, col=1)
    # yrange for H response (1,2) is equal to yrange for M-H (1,1)
    if max_Hresponse > max_Mresponse:
        max_Hresponse = max_Mresponse        
    fig.update_yaxes(range=[-max_Hresponse, max_Hresponse], row=1, col=2)
    fig.update_yaxes(range=[-max_Hresponse, max_Hresponse], row=1, col=1)

    return fig

def generate_langevin_data(H_max=100e-3/mu_0):
    """
    Generate Langevin function data for plotting.
    """
    H_L = np.linspace(-1*H_max, H_max, 500)
    M_L = langevin(H_L)
    return H_L, M_L

def Response_several_waveforms(Input_waveform="sinusoidal"):
    """
    Main function to generate and plot Langevin function response.
    """
    # Parameters
    omega_drive = 2 * np.pi * 1e3  # 1kHz
    Hamp_drive = 20e-3 / mu_0  # 1mT/mu_0
    Hdc_drive = 0e-3 / mu_0  # 0mT/mu_0
    T_cycles = 3
    Npts_cycles = 500

    # Generate applied field based on the selected waveform
    if Input_waveform == "sinusoidal":
        t, Hac, Fs = generate_applied_field(omega_drive, Hamp_drive, Hdc_drive, T_cycles, Npts_cycles)
    elif Input_waveform == "square":
        t, Hac, Fs = generate_square_wave(1e3, 20e-3 / mu_0, 0, T_cycles, Npts_cycles)
    elif Input_waveform == "triangle":
        t, Hac, Fs = generate_triangle_wave(1e3, 20e-3 / mu_0, 0, T_cycles, Npts_cycles)
    elif Input_waveform == "chirp":
        t, Hac, Fs = generate_chirp_signal(1e3, 2e3, T_cycles * 2 * np.pi / omega_drive, Npts_cycles, amplitude=20e-3 / mu_0, offset=0, method='linear')
    elif Input_waveform == "morlet":
        t, Hac, Fs = generate_morlet_wavelet(1e3, 0.1)
    elif Input_waveform == "mexican_hat":
        t, Hac, Fs = generate_mexican_hat_wavelet(0.1)
    elif Input_waveform == "haar":
        t, Hac, Fs = generate_haar_wavelet()
    else:
        raise ValueError("Invalid waveform type")

    # Generate Langevin function response
    H_resp = generate_response(Hac)

    # Generate Langevin function data
    H_L, M_L = generate_langevin_data(np.max(Hac)*1.2)

    # FFT processing
    f_H, A_H, _ = compute_fft(Hac, Fs)
    f_M, A_M, _ = compute_fft(H_resp, Fs)

    # Create figure
    fig = create_figure(t, Hac, H_resp, H_L, M_L, f_H, A_H, f_M, A_M)
    
    return fig

2025-05-15

【技術メモ】PythonでMPI Physics Demo を再現してみた-Part 1

MIT より Educational Simulation として MPI Physics Demo が公開されていた。

github.comMATLABで作成されており、アニメーションで視覚的に原理の理解が出来る。

汎用性を高めるために、Pythonで同シミュレーションの再現に取り組むことを考えている。アニメーションまでは出来ていないが、基本原理の部分の再現コードを作成した。

実行結果を以下のプロットに示す。
第二象限は粒子のＭ-H曲線
第三象限は印加磁場波形（横軸が磁場強度、縦軸を時間としている）
第一象限は粒子磁化応答波形（横軸が時間、縦軸は磁場強度）
第四象限は印加磁場、粒子磁荷応答のFFT結果（最大値で正規化した結果）

コードは下記（清書していないため、お見苦しい点が多々あると思われる）。

import numpy as np
from dash import Dash, html, dcc
import plotly.graph_objs as go

##################################################################
################# Langevin function ##############################
##################################################################
def langevin(x, c=1, Ms=5e5, D=20e-9, kb=1.380649e-23, T=300): 
    # Langevin function
    # x: strength of the magnetic field 'H'
    # kb: Boltzmann constant 1.380649e-23  # J/K
    # Ms: saturation magnetization
    # D: diameter of the particle
    # c: concentration
    # T: temperature
    mu_0 = 4*np.pi*1e-7  # permeability of free space
    V = np.pi * D**3 / 6  # volume of the particle
    m = Ms * V  # magnetic moment
    # Langevin parameter
    beta = mu_0 * m / (kb * T)  # Langevin parameter
    coefficient = c*m
    bH = beta * x
    """Langevin function L(x) = coth(x) - 1/x"""
    y = np.zeros_like(bH)
    for i in range(len(bH)):
        if bH[i] == 0:
            y[i] = 0
        else:
            y[i] = coefficient * (1/np.tanh(bH[i]) - 1/bH[i])
    return y

######################################################################################
######################### Langevin function response ####################################
######################################################################################
# Generate Applied Field data
mu_0 = 4*np.pi*1e-7  # permeability of free space
omega_drive = 2*np.pi*1e3# 1kHz
Hamp_drive = 20e-3 / mu_0 # 1mT/mu_0
Hdc_drive = 0e-3 /mu_0 # 0mT/mu_0

# Generate Sweep data
T_cycles = 3
Npts_cycles = 500
Fs = 1/(2*np.pi/(omega_drive)*T_cycles/Npts_cycles)
t = np.linspace(0, 2*np.pi/(omega_drive)*T_cycles, Npts_cycles) # Time vector 3 cycles
Hac_drive = Hamp_drive * np.sin(omega_drive * t) + Hdc_drive

# Generate Langevin function response
H_response = langevin(Hac_drive)

# Generate data for Langevin function
H_Langevin = np.linspace(-100e-3, 100e-3, 500) / mu_0 # Avoid x=0 to prevent division by zero
M_Langevin = langevin(H_Langevin)

######################################################################################
####### FFT Result of Applied Field data and Langevin function response ##############
######################################################################################
# FFT of Hac_drive (correlation of amplitude by muliplying by 1/Number of points at DC component or 2/Number of points at otherwise)
Hac_drive_fft = np.fft.fft(Hac_drive) * 2 / len(Hac_drive)
Hac_drive_freq = np.fft.fftfreq(len(Hac_drive), d=1/Fs)
Hac_drive_fftamplitude = np.abs(Hac_drive_fft) * 2 / len(Hac_drive)
Hac_drive_fftamplitude[0] = Hac_drive_fftamplitude[0] / 2 # DC component
Hac_drive_fftphase = np.angle(Hac_drive_fft)

# FFT of H_response (correlation of amplitude by muliplying by 2/Number of points)
H_response_fft = np.fft.fft(H_response) * 2 / len(H_response)
H_response_freq = np.fft.fftfreq(len(H_response), d=1/Fs)
H_response_fftamplitude = np.abs(H_response_fft) * 2 / len(H_response)
H_response_fftamplitude[0] = H_response_fftamplitude[0] / 2 # DC component
H_response_fftphase = np.angle(H_response_fft)



######################################################################################
######################### Create figure ##############################################
######################################################################################
fig = make_subplots(
    rows=2, cols=2,
    subplot_titles=('M(H)', 'H_response(t)', 'H_applied(t)'),
    vertical_spacing=0.1
)
# Plot Langevin function
fig.add_trace(go.Scatter(x=H_Langevin*mu_0, y=M_Langevin, mode='lines', name='M(H)', line=dict(color='black', width=1, dash='dash')), row=1, col=1) # dashed line
fig.add_trace(go.Scatter(x=Hac_drive*mu_0, y=H_response, mode='lines', name='Sweep', line=dict(color='green', width=2)), row=1, col=1)
# FFT of Hac drive
fig.add_trace(go.Scatter(y=t, x=Hac_drive*mu_0, mode='lines', name='H(t)', line=dict(color='blue')), row=2, col=1)
# FFT of H_response
fig.add_trace(go.Scatter(x=t, y=H_response*mu_0, mode='lines', name='M(t)', line=dict(color='red')), row=1, col=2)
# FFT of H_response and Hac_drive
# fig.add_trace(go.Scatter(x=H_response_freq[0:Npts_cycles//2], y=H_response_fftamplitude[0:Npts_cycles//2]*mu_0, mode='lines', name='FFT(H_response)', line=dict(color='red')), row=2, col=2)
# fig.add_trace(go.Scatter(x=Hac_drive_freq[0:Npts_cycles//2], y=Hac_drive_fftamplitude[0:Npts_cycles//2]*mu_0, mode='lines', name='FFT(H_applied)', line=dict(color='blue')), row=2, col=2)
fig.add_trace(go.Scatter(x=H_response_freq[0:Npts_cycles//2], y=H_response_fftamplitude[0:Npts_cycles//2]/np.max(H_response_fftamplitude[0:Npts_cycles//2]), mode='lines', name='FFT(H_response)', line=dict(color='red')), row=2, col=2)
fig.add_trace(go.Scatter(x=Hac_drive_freq[0:Npts_cycles//2], y=Hac_drive_fftamplitude[0:Npts_cycles//2]/np.max(Hac_drive_fftamplitude[0:Npts_cycles//2]), mode='lines', name='FFT(H_applied)', line=dict(color='blue')), row=2, col=2)

# For row 1 and col 1, x-label is 'H/mu_0 [T]', y-label is 'M(H)/mu_0 [T]'
# For row 2 and col 1, y-label is 't [s]', x-label is 'H_applied/mu_0 [T]'
# For row 1 and col 2, x-label is 't [s]', y-label is 'H_response/mu_0 [T]'
fig.update_xaxes(title_text='H [T]', row=1, col=1)
fig.update_yaxes(title_text='M(H) [A/m]', row=1, col=1)
fig.update_yaxes(title_text='t [s]', row=2, col=1)
fig.update_xaxes(title_text='H_applied [T]', row=2, col=1)
fig.update_xaxes(title_text='t [s]', row=1, col=2)
fig.update_yaxes(title_text='M(t) [A/m]', row=1, col=2)
fig.update_xaxes(title_text='Frequency [Hz]', row=2, col=2)
fig.update_yaxes(title_text='FFT [A/m]', row=2, col=2)
# For row 2 and col 2, xrange is [0, Fs/2]
fig.update_xaxes(range=[0, np.log10(Fs/2)], row=2, col=2)
fig.update_xaxes(type='log', row=2, col=2)  # log scale for x-axis in subplot (2,2)
fig.update_yaxes(type='log', row=2, col=2)  # log scale for y-axis in subplot (2,2)

fig.update_layout(
    height=800, width=1200,
    title_text='Langevin Function Response',
    showlegend=False,
    margin=dict(l=40, r=40, t=60, b=40)
)
# Create Dash app
app = Dash(__name__)
app.layout = html.Div([
    html.H2("Langevin Function Response"),
    dcc.Graph(
        id='subplot-graph',
        figure=fig
    )
])

if __name__ == "__main__":
    app.run_server(debug=True)

2025-05-04

PythonでExcel/CSVデータのFFT処理自動化

技術メモ

はじめに

実験やパラメータスタディで得られた時系列データを大量にFFT（高速フーリエ変換）処理したい場面は、研究・開発の現場では珍しくない。実験結果の時系列データを複数種類・複数条件でFFT処理しなければならない場面、パラメータスタディをした解析結果の複数データをFFT処理しなければならない場面、...　etc. など、筆者もこれまで多くの場面で「大量のデータをFFT処理」することが要求される場面に出くわしてきた。
しかし、エクセルのアドインだと処理が重たいし、MATLABやPythonで処理するとなるとfftするところまではよいが、振幅の係数処理をどうするだの、fft結果をexcelファイル等に出力する処理をどうするだの、何かと面倒なことが多かった。
今回、そこで本稿では、CSV／ExcelファイルのデータをPythonで一括FFTし、振幅補正や窓関数適用済みの結果をそのままCSV／Excelに出力するプログラムを作成した。高度なプログラミング知識がなくても、Python環境さえあればすぐに実行できるよう設計してある。

本プログラムのポイント

直感的な対象データ範囲の指定が可能
ファイルの選択はファイルダイアログにてGUI上で指定。
エクセルのセル範囲（D2:F18 など）で対象データ範囲指定を行う。
振幅補正済みの結果を出力
計算結果の振幅値を1/N補正する処理や、窓関数を適用した際の振幅補正処理済みの結果を出力する。
出力形式、方式を複数のバリエーションから選択できる。
振幅形式（「振幅/位相形式」または「cos振幅/sin振幅形式」）の選択や、
データの出力順序（振幅データを全て出力したのちに位相データを出力、や、元データごとに振幅、位相を交互に出力）の方式選択を行うことが出来る
Jupyter Notebook で処理コードの関数を先に実行し、実行コードを後のブロックで実行するような使い方をすると便利。
処理コードは関数化しており、実行コードは1行で済むように実装している。先に関数をロードしておいて、後のブロックで実行コードを実行する使い方が出来る。

引数仕様と使い方例

〇実行コード（引数値は一例）

fft_excel_with_dialog(
    dt=0.01,            # サンプリング周期（秒）
    start_row=2,        # データ開始行
    start_col='B',      # データ開始列
    end_col='E',        # データ終了列
    output_format=1,    # 出力形式（省略可、デフォルトは1）
    sheet_name='Sheet1',# Excelシート名（省略可）
    window_name='hanning'# 窓関数の種類（省略可、デフォルトは矩形窓）
)

引数名	説明	例
`dt`	サンプリング周期（秒）。等間隔サンプリング前提。例：1 MHzデータなら `1.0e-6` を指定する。	`0.01`
`start_row`	FFT対象データの開始行番号（1-index）。ヘッダ行がある場合はヘッダ分だけずらす。	`7`
`start_col`	FFT対象データの開始列（A–Z、複数文字可）。	`'D'`
`end_col`	FFT対象データの終了列（A–Z、複数文字可）。	`'I'`
`output_format`	1–4 の出力形式。省略時は `1`（振幅→位相の順）。詳細は次節参照。	`2`
`sheet_name`	複数シートを含むExcelファイルを指定する場合のシート名。省略時は先頭シートを利用。	`'waveform2'`
`window_name`	窓関数名。省略時は矩形窓（`rectangular`）となり、窓処理なしと同等。`hanning`、`hamming`、`blackmanharris` 等を利用可能。

範囲指定について

start_row, start_col, end_col, sheet_nameではFFT対象データの範囲、シート名を指定する。例えば、以下の 'waveform2' シートの範囲（D7:I190）のデータを処理したい場合は次のように指定する。

start_row=7, start_col='D', end_col='I', sheet_name='waveform2'

`output_format` の4パターン

出力形式・順序を指定する。省略した場合は、振幅/位相形式で、「振幅、振幅、…、位相、位相、…、位相」という順序で出力される。

引数値	出力形式	出力順序	例
1	振幅/位相	全データの振幅を出力したのちに、全データの位相を出力	データ1の振幅、データ2の振幅、... 、データMの振幅、データ1の位相、データ2の位相、... 、データMの位相
2	振幅/位相	振幅、位相の順番でデータごとに出力	データ1の振幅、データ1の位相、データ2の振幅、データ2の位相、... 、データMの振幅、データMの位相
3	cos振幅/sin振幅	全データの振幅を出力したのちに、全データの位相を出力	データ1のcos振幅、データ2のcos振幅、... 、データMのcos振幅、データ1のsin振幅、データ2のsin振幅、... 、データMのsin振幅
4	cos振幅/sin振幅	cos振幅、sin振幅の順番でデータごとに出力	データ1のcos振幅、データ1のsin振幅、データ2のcos振幅、データ2のsin振幅、... 、データMのcos振幅、データMのsin振幅

窓関数の選び方

FFT前に窓関数を適用すると、スペクトル漏れ（リーク）を抑制できる。代表的な例は以下の通り：

矩形窓 (rectangular)：窓処理なしと同等
ハニング窓 (hanning)：バランスの良い選択
ハミング窓 (hamming)：主ローブ幅を最小化
ブラックマン–ハリス窓 (blackmanharris)：サイドローブ抑制重視
フラットトップ窓 (flattop)：振幅精度を最重視

今回、'hanning', 'hamming', 'blackmanharris', 'flattop' など有名どころの窓関数は選択できるように実装している。（コード中に選択可能な窓関数が列挙されている。）
下記に記載のように、適用対象の性質に応じて使い分けるとよいと考えられる。

実行コードおよび処理コード

実行コードおよび処理部のコードは以下に記載した。
実行コードは1行で完結しており、引数値のみ対象データに合わせて変更し、実行すれば処理が開始される。

Jupyter Notebook の先頭ブロックで処理コードを先に実行しておき、実行コードを別のコマンドブロックで実行するような使い方をすると便利かもしれない。

〇実行コード（引数値は一例）

fft_excel_with_dialog(dt=0.01,start_row=2,start_col='B',end_col='E',output_format=1,sheet_name='Sheet1',window_name='hanning')

〇処理コード
主な処理フローおよび、コード詳細を次の通りである。

処理コードの主な処理フロー（概要）

ファイル選択ダイアログを起動
引数チェック
CSV／Excel読み込み（pandas＋openpyxl）
データ抽出（開始行・開始列～終了列）
窓関数生成および補正係数計算
FFT実行
- 実数部／虚数部を正規化し、cos成分・sin成分を分離
- 振幅・位相を算出
結果DataFrameの組み立て
CSV保存またはExcelへの追加出力

import tkinter as tk
from tkinter import filedialog
import os
import numpy as np
import pandas as pd
from openpyxl import load_workbook
from openpyxl.utils import get_column_letter, column_index_from_string


def fft_excel_with_dialog(
    dt: float,
    start_col: str,
    start_row: int,
    end_col: str,
    output_format: int = 1,
    sheet_name: str = None,
    window_name: str = 'rectangular'
):
    """
    Excelの任意範囲を選択してFFTを行い、結果を新規シートに追加

    Args:
      dt: サンプリング間隔（秒）
      start_row: FFT対象データの開始行番号 (1-index)
      start_col: FFT対象データの開始列 (例: 'F','BG')
      end_col: FFT対象データの終了列 (例: 'J','BH')
      output_format: 1-4 の出力形式
        1: all amps, then all phases
        2: amp,phase per data interleaved
        3: all cos amps, then all sin amps
        4: cos,sin per data interleaved
      sheet_name: Excelシート名（省略可）
      window_name: 'rectangular','hanning','hamming',...        
    """
    # 1) ファイル選択ダイアログ
    root = tk.Tk()
    root.attributes('-topmost', True)
    root.withdraw()
    path = filedialog.askopenfilename(
        title="FFTを行うExcelファイルを選択",
        filetypes=[("CSV/Excel files","*.csv;*.xlsx;*.xls;*.xlsm")
        ]
    )
    root.destroy()
    if not path:
        print("ファイルが選択されませんでした。処理を終了します。")
        return

    # 2) データ読み込み前の準備
    # 2-1) 引数の検証
    if dt <= 0:
        raise ValueError("dt must be > 0")
    if start_row < 1:
        raise ValueError("start_row must be >= 1")
    if start_col < 'A' or start_col > 'Z':
        raise ValueError("start_col must be a valid column letter (A-Z)")
    if end_col < 'A' or end_col > 'Z':
        raise ValueError("end_col must be a valid column letter (A-Z)")
    if start_col >= end_col:
        raise ValueError("start_col must be < end_col")
    if output_format not in [1, 2, 3, 4]:
        raise ValueError("output_format must be 1-4")
    
    # 2-2) データ読み込み前の準備
    # ヘッダ行を指定してデータを読み込む
    if start_row == 1:
        header_row = None
    elif start_row < 1:
        raise ValueError("start_row must be >= 1")  
    else:
        header_row = start_row - 2
    # 拡張子取得
    base, ext = os.path.splitext(path)
    ext = ext.lower()    

    # 3) データ読み込み
    # データ読み込み
    if ext == '.csv':
        df = pd.read_csv(path, header=header_row)
        sheet_in = None
    else:
        # Excel読み込み
        book = load_workbook(path, read_only=True)
        if sheet_name is not None:
            # 指定されたシート名が存在するか確認
            if sheet_name not in book.sheetnames:
                raise ValueError(f"Sheet '{sheet_name}' does not exist in the Excel file.")
            sheet_in = sheet_name
        else:
            # 最初のシートを使用
            sheet_in = book.sheetnames[0]
        df = pd.read_excel(path, sheet_name=sheet_in, engine='openpyxl', header=header_row)

    # データの行数と列数を取得    
    col_start = column_index_from_string(start_col)
    col_end = column_index_from_string(end_col)
    data = df.iloc[start_row-1:, col_start-1:col_end]
    if header_row is None:
        data.columns = [get_column_letter(i) for i in range(col_start, col_end+1)]

    # 4) FFT準備
    N = len(data)
    freq = np.fft.fftfreq(N, d=dt)

    # 窓関数選択
    # ─ 窓関数の選択 ───────────────────────────────────────
    if window_name == 'blackmanharris':
        window = np.blackmanharris(N)  # ブラックマン-ハリス窓
    elif window_name == 'bartlett':
        window = np.bartlett(N)  # バートレット窓
    elif window_name == 'kaiser':
        beta = 14.0  # ベータ値（カイザー窓の形状を決定するパラメータ）
        window = np.kaiser(N, beta)  # カイザー窓
    elif window_name == 'flattop':
        window = np.flattop(N)  # フラットトップ窓
    elif window_name == 'blackmanharris4':
        window = np.blackmanharris(N, sym=False)  # ブラックマン-ハリス4窓
    elif window_name == 'nuttall':
        window = np.nuttall(N)  # ナッタール窓
    elif window_name == 'tukey':
        alpha = 0.5  # タキ窓の形状を決定するパラメータ
        window = np.tukey(N, alpha)  # タキ窓
    elif window_name == 'cosinebell':
        window = np.cos(np.linspace(0, np.pi, N))**2  # コサインベル窓
    elif window_name == 'triangular':
        window = np.triang(N)  # 三角窓
    elif window_name == 'parzen':
        window = np.parzen(N)  # パルゼン窓
    elif window_name == 'hanning':
        window = np.hanning(N)  # ハニング窓
    elif window_name == 'hamming':
        window = np.hamming(N)  # ハミング窓
    elif window_name == 'rectangular':
        window = np.ones(N)
    else:
        raise ValueError(f"Unknown window name: {window_name}")
    # 窓関数適用に対する振幅補正値
    window_correctionfactor = N / np.sum(window)  # 窓関数の積分値で割る    


    # 5) FFT実行
    cols = data.columns.tolist()
    results = {}
    for col in cols:
        y_win = data[col].astype(float).to_numpy() * window
        Y = np.fft.fft(y_win)
        # cosine成分とsine成分を分離
        Re = Y.real
        Im = Y.imag
        re_norm = np.empty_like(Re)
        im_norm = np.empty_like(Im)
        re_norm[0] = Re[0] / N
        re_norm[1:] = Re[1:] * 2 / N
        im_norm[0] = Im[0] / N
        im_norm[1:] = Im[1:] * 2 / N
        # 振幅と位相を計算
        A = np.sqrt(re_norm**2 + im_norm**2)
        phase = np.angle(Y, deg=True)
        # 結果を辞書に格納
        results[col] = {
            'amp': window_correctionfactor * A,
            'phase': phase,
            'cos_amp': window_correctionfactor * re_norm,
            'sin_amp': window_correctionfactor * im_norm
        }

    # 6) 出力DataFrame作成
    out = pd.DataFrame({'Frequency': freq})
    if output_format == 1:
        for col in cols:
            out[f'{col}_amp'] = results[col]['amp']
        for col in cols:
            out[f'{col}_phase'] = results[col]['phase']
    elif output_format == 2:
        for col in cols:
            out[f'{col}_amp'] = results[col]['amp']
            out[f'{col}_phase'] = results[col]['phase']
    elif output_format == 3:
        for col in cols:
            out[f'{col}_cos_amp'] = results[col]['cos_amp']
        for col in cols:
            out[f'{col}_sin_amp'] = results[col]['sin_amp']
    elif output_format == 4:
        for col in cols:
            out[f'{col}_cos_amp'] = results[col]['cos_amp']
            out[f'{col}_sin_amp'] = results[col]['sin_amp']
    else:
        raise ValueError("output_format must be 1-4")

    
    # 7) 出力ファイル名の決定
    out_path = f"{base}_FFT{ext}"    

    # 8) 出力ファイルの保存
    if ext == '.csv':
        # 8-1) CSVファイルの場合
        out_path = f"{base}_FFT.csv"
        out.to_csv(out_path, index=False)
        print(f"FFT結果をCSV保存: {out_path}")
    else:
        # 8-2) Excelファイルの場合
        # シート名重複回避
        base_name = 'FFT'
        out_sheet = base_name
        suffix = 1
        while out_sheet in book.sheetnames:
            out_sheet = f"{base_name}_{suffix}"
            suffix += 1        
        # 新規シートに出力
        with pd.ExcelWriter(out_path, engine='openpyxl', mode='a', if_sheet_exists='new') as writer:
            out.to_excel(writer, sheet_name=out_sheet, index=False)

        print(f"FFT結果をシート'{out_sheet}'に追加しました: {path}")

おわりに

本プログラムは、大量FFT処理後のわずらわしい後処理を自動化し、作業効率の向上に寄与することを目的としており、皆様の作業効率化の一助となれば幸いである。バグや機能追加の要望、ブラッシュアップした改造案（コード）などあれば、遠慮なくコメントを頂けると嬉しい限りである。

2025-05-03

【論メ】Low-Concentration Magnetic Particle Spectroscopy Using Gradiometric Receive Coil-Coupled Magnetoresistive Sensor

論文

〇論文情報：

・タイトル

　Low-Concentration Magnetic Particle Spectroscopy Using Gradiometric Receive Coil-Coupled Magnetoresistive Sensor

・著者

　S. B. Trisnanto 他

　横浜国立大学、TDKのチーム

・会議、ジャーナル

　IEEE TRANSACTIONS ON MAGNETICS, VOL. 59, NO. 11 (2023)

〇アブストラクトの日本語訳

MPIの臨床応用を実現するには、磁気ナノトレーサーを安全に投与できる範囲内で検出できるように、スキャナーを高感度に設計する必要がある。ゼロ次元MPIスキャナーとして、強度の大きい正弦波磁場を印加することでトレーサーの複雑な挙動を特徴づけるためにMPSを使用することが多い。しかし、トレーサー濃度が2µgFe/mLを下回ると、磁化信号がとても小さくなり、細胞外基質の電磁現象に支配されてしまう。今回、我々は典型的なグラジオメーター受信コイルに対してMRセンサーを結合させることによって、10kHz・25mT/µ0 励磁場印加時にMPSの感度を向上させることが出来ることを実証した。この磁気測定系では、70ngFeを含んだ0.1mL希釈Resovistサンプルの信号を弁別することに成功した。また、本研究の試験では、磁場に誘起されて水をベースとした溶媒から生じたバックグラウンド信号が非常に大きいことを確認した。イオン濃度に応じて、極性を有する溶液からは反磁性および、渦電流による影響の両方が現れることが分かった。たとえ高調波成分が理論上消失したとしても、低濃度サンプルの磁化信号は基底周波数において、媒体（背景溶液）の信号と比べて位相や振幅に差異を示していた。その高い感度を活かせば、MRセンサーを用いたMPSシステムは、MPI技術に加えて液相バイオセンシングへの応用も期待できる。

〇課題、目的と何をどうやって検証したのか？

・課題

低濃度トレーサーを用いたMPSを生体センシングへの適用可能性を検討している。

MPSのSNR向上のための従来手法では、グラジオコイルの感度向上や信号ミキシングによるノイズ低減が用いられるが、

磁気ヒステリシスの再構成が困難であり、磁気粒子の緩和挙動を測定することが出来なかった。

・目的

SNR向上のためにTMRセンサを用いたグラジオコイルの広帯域ノイズを補正する装置案の実証と、

高SNR測定系を用いた磁気ナノ粒子の電磁気特性評価

①SNRの改善効果の実証

②溶媒の電磁現象が磁気ヒステリシスに与える影響の評価

③粒子濃度と高次信号強度のリニアリティ評価

・検証

①グラジオの検出信号に伴う誘導電流を二次コイルに流し、二次コイル磁束をTMRセンサで検出する構成を実装。

　7ugFe, 140ngFe, 水, B.G信号をTMRセンサ有無で比較。

　波形のノイズが除去された様子および、ノイズ成分の振幅が低下していることを確認することで実証。

　ノイズが低減されたため、反磁性効果により水において励磁場よりも位相遅れが生じることも確認。

②.1：低濃度サンプルでは高次信号が低減し、磁化ヒステリシス（M-H）が線形的な挙動となった。

　測定電圧の積分からMを計算。70ugFe~70ngFe/0.1mL Risovistサンプルを測定し、M-H曲線から確認。

②.2：イオン化溶媒の渦電流により信号強度の増加と位相進みを確認。水の反磁性による励磁場との位相反転応答を確認。

　水、様々なモル濃度（0.07~6.00M）で導電率を調整したNaCl溶液、導電率が非常に低いシリコンオイルの基本波信号応答

　測定結果を比較。測定電圧の積分からMを計算。70ugFe~70ngFe/0.1mL Risovistサンプルを測定し、M-H曲線から確認。

③70ugFe~70ngFe/0.1mL Risovistサンプルを測定し、1,3,5,…,61次信号と濃度の関係を確認。

　3次,5次,7次までは最低濃度(70ngFe)まで濃度とリニアリティがあることを確認。

〇何を新規性として強調しているか

・TMRセンサを用いてグラジオコイルの広帯域ノイズを低減した測定を実現し、受信系のSNR向上を達成

・25mT/µ0, 10kHzの励磁系を用いて、7µgFe/mL 以下の濃度のRisovistサンプルからの信号検出を達成

・基本波信号を分析することで、バックグラウンド信号から溶媒の反磁性影響、渦電流影響を特定できる

・鉄濃度が非常に低くなると、磁気ヒステリシスが線形挙動を示すことを実験的に明らかにした

〇技術、手法、アイディアなどでキーとなっている点、凄いと思った点

・システム全体

・ソレノイドコイル（Tx）の内側にグラジオコイル（Rx）を配置。

　グラジオコイルはボルテージフォロワー回路（アンプ回路）を介して磁気シールドチャンバー内の二次コイルに接続。

　二次コイルはグラジオコイル検出信号に伴う誘導電流に比例した磁束を生じ、TMRセンサがそれを検出する。

　TMRセンサの検出信号が最終的なアウトプット信号。

・励磁コイル

　300ターン、長さ10mm、10.2mT/A、25mT/µ0 @ 10kHz

・グラジオコイル

　50ターン、外径Φ10mm

・二次コイル

　300ターン、Φ3mm

・TMRセンサ

　TDK Nivio xMRセンサ、12x12x74mm、87µV/nT@5V駆動時

　1MHzカットオフのデジタルLPFを実装

・磁気シールドチャンバー

　2層のPyチューブ＆Pyシートで構成、外部磁場を0.7nTまで抑制

・グラジオコイルの検出信号を再度磁束に変換し、TMRセンサ検出を行うことでグラジオ起因の広帯域ノイズを消去していた。

・イオン性溶媒、反磁性溶媒の基本波位相に着目し、比較実験から反磁性影響、渦電流影響を特定していた。

・測定電圧信号の積分＆正規化から磁化特性を評価していた。

〇この論文の限界

・SNRを定量評価していない。

・BG信号の傾向の周波数依存性、励磁強度依存性で評価していない。

〇Discussionでの興味深い仮説や解釈

・イオン化溶媒中の磁性粒子を測定する場合は、バックグラウンド信号のイオン濃度依存性を考慮する必要があること

反磁性溶媒（シリコンオイル）では励磁場に対して磁化応答の位相が反転＆強度が若干低くなる

イオン化溶媒（NaCl溶液）ではモル濃度が高いと渦電流影響が大きくなり、励磁場と磁化応答の位相差が小＆強度が若干高まる

・低濃度では、Langevin関数によるモデリングで説明できるため、濃度と高次信号強度が比例する解釈

Langevin関数による超常磁性粒子のモデル化は磁気ポテンシャル（µ0m）と熱擾乱（kbT)のバランスを現象論的に説明したもの磁気的な相互干渉や、磁気異方性の影響を考慮していない式とのこと。測定結果のm3, m5, m7, m9は濃度と信号強度が比例しているため、この解釈があっている

〇あなたの気づき、自分への活用案

・センサのドリフト影響はあるのか？

〇次に読むべき論文

・参考文献[12]：グラジオコイルの詳細が記載されている。

・参考文献[11]：今回と同様の原理の磁場計測事例が紹介されている。

・参考文献[19]：相互干渉や異方性の影響を考慮した超常磁性粒子の特性が記述されているかもしれない。

・参考文献[8]~[10]：その他のMPS構成を学べる可能性あり。

〇その他参考情報

★ボルテージフォロワー

　利得1で、入力電圧をそのまま出力する。

　オペアンプを介して、高い入力インピーダンスと低い出力インピーダンスを確保できるため、バッファとして用いられる。

引用元：https://toragi.cqpub.co.jp/wp-content/uploads/p092-19.pdf

★出力インピーダンス、入力インピーダンスについて

　「出力インピーダンス小」＆「入力インピーダンス大」の方が電圧損失を抑えて信号伝送ができる！

https://detail-infomation.com/output-impedance-input-impedance/#google_vignette

とある大学院生の趣味備忘録

Python, マイコンいじり, 日々の呟きなど

Plotlyでエクセル上のデータをプロットするスクリプトを作成した。

はじめに

Excel のグラフ描画機能における課題

1. プロットデータの登録方法が非効率的

2. デフォルトテンプレートの設定が変更できない

3. 描画範囲の調整に手間がかかる

Python＋Plotly で散布図を自動生成するスクリプト

使い方例1

使い方例②：追加データを重ねてプロット

関数スクリプト

Pythonでタスク分類ツールを作った

タスク分類ツール

機能

メトリクス

ジャンル

使用方法

開発環境

インストール

実行方法

ノッチフィルタ（バンドストップフィルタ）の勉強備忘録

【技術メモ】PythonでMPI Physics Demo を再現してみた-Part 2

【技術メモ】PythonでMPI Physics Demo を再現してみた-Part 1

PythonでExcel/CSVデータのFFT処理自動化

はじめに

本プログラムのポイント

引数仕様と使い方例

範囲指定について

`output_format` の4パターン

窓関数の選び方

実行コードおよび処理コード

処理コードの主な処理フロー（概要）

おわりに

【論メ】Low-Concentration Magnetic Particle Spectroscopy Using Gradiometric Receive Coil-Coupled Magnetoresistive Sensor

はじめに

Excel のグラフ描画機能における課題

1. プロットデータの登録方法が非効率的

2. デフォルトテンプレートの設定が変更できない

3. 描画範囲の調整に手間がかかる

Python＋Plotly で散布図を自動生成するスクリプト

使い方例1

使い方例②：追加データを重ねてプロット

関数スクリプト

タスク分類ツール

機能

メトリクス

ジャンル

使用方法

開発環境

インストール

実行方法

はじめに

本プログラムのポイント

引数仕様と使い方例

範囲指定について

output_format の4パターン

窓関数の選び方

実行コードおよび処理コード

処理コードの主な処理フロー（概要）

おわりに

`output_format` の4パターン