S transform

A: 良い質問です。 ウェーブレット変換（Wavelet Transform, WT）の「対抗馬」として提案された S変換（Stockwell Transform, ST） は、時間-周波数解析における重要な手法の一つです。 以下では、S変換の 原理・数理的特徴・ウェーブレット変換との比較・応用例 の順で体系的に説明します。

---

## 🧭 1. 背景と目的

フーリエ変換（FT）は信号を周波数成分に分解しますが、/時間情報が失われる/ため、非定常信号（時間とともに周波数が変化する信号）の解析には不向きです。 そこで登場したのが以下の手法です：

手法	時間・周波数の扱い	主な特徴
----------------	--------------------	----------------
短時間フーリエ変換 (STFT)	固定ウィンドウ	解像度が固定、柔軟性がない
ウェーブレット変換 (WT)	可変ウィンドウ	周波数に応じたスケーリングが可能
S変換 (ST)	フーリエ変換＋周波数依存ガウス窓	FTとWTの両方の長所を併せ持つ

S変換は、/STFTの窓長を周波数に依存して変化させる/ことにより、 低周波成分には広い窓、高周波成分には狭い窓を適用します。 これにより、/時間分解能と周波数分解能のバランスを自動的に調整/できます。

---

## ⚙️ 2. 数理的定義と原理

S変換は、信号 ( x(t) ) に対して以下のように定義されます：

[ S(τ, f) = ∫_-∞^∞ x(t) , \frac{|f|}{\sqrt{2\pi}} e^{-\frac{(t - \tau)^2 f^2}{2}} e^{-i 2 π f t} , dt ]

### 各要素の意味：

• ( τ )：時間位置（タイムシフト）
• ( f )：周波数
• ( e^{-\frac{(t - \tau)^2 f^2}{2}} )：ガウス窓（周波数依存）
• ( e^{-i 2π f t} )：フーリエカーネル

つまり、

• /フーリエ変換のように全周波数を扱いつつ/、
• /ガウス窓の幅を周波数に応じて変える/ことで、 高周波ほど時間分解能が上がり、低周波ほど周波数分解能が上がる。

この構造から、S変換は 「/ウェーブレット変換のようなスケーリング性/を持つが、/フーリエ変換のように位相を保つ/」 という特徴を持ちます。

---

## 🧩 3. ウェーブレット変換との比較

特徴	ウェーブレット変換 (WT)	S変換 (ST)
-------	------------------	-------------------
母関数	任意（Morlet, Haarなど）	ガウス型（Morletに近い）
解析基底	スケーリング・シフト	周波数依存ガウス窓＋FT
位相情報	通常は失われる	保持される（複素値）
再構成の容易さ	変換の種類による	容易（逆S変換が明示的）
計算負荷	やや低い	やや高い（FT利用）
適用信号	非定常信号、急変	地震波、音声、医用信号など非定常データ
特徴	スケール解析に優れる	時間-周波数解析＋位相保持

特に、/位相情報を保持できる点/がS変換の強みです。 ウェーブレット変換はしばしば振幅解析に強い一方で、S変換は/波形の位相進行/や/群遅延/の解析に適しています。

---

## 🧮 4. 計算面での工夫

S変換は定義のままでは積分計算が重いため、 実際には/フーリエ領域での高速化（FFTベース）/がよく用いられます。

離散信号 ( x[n] ) に対して：

1. ( X[k] = \text{FFT}(x[n]) )
2. 各周波数 ( f_k ) ごとにガウス窓を適用
3. 逆FFTを取って時間領域に戻す

この方法を「/高速S変換（Fast S-transform）/」と呼びます。

---

## 🌍 5. 応用例

S変換は、/時間変化のある周波数構造を持つ信号の解析/に特に有効です。

### 代表的な応用領域：

分野	用途
--------	-----------------------
地球物理・地震学	地震波の瞬時周波数解析、地層の時間-周波数構造
医用信号解析	EEG・ECGのイベント検出、異常波形識別
音響解析	瞬時スペクトルの追跡、音の特徴抽出
電力システム	瞬時電力品質の解析、異常検出
機械振動診断	ベアリングの異常検知、変動周波数振動の追跡

---

## 🔍 6. まとめ（要点整理）

項目	内容
---	-----------------------------------------------
名称	S変換（Stockwell Transform）
提案者	R. G. Stockwell, L. Mansinha, R. P. Lowe (1996)
本質	フーリエ変換に周波数依存ガウス窓を導入した時間-周波数変換
特徴	フェーズ保持・周波数依存ウィンドウ・再構成容易
強み	非定常信号の解析、瞬時位相・瞬時周波数の抽出
弱み	計算コストがやや高い（FFTで軽減可能）

import numpy as np
import matplotlib.pyplot as plt

def fast_s_transform(x, fs):
    """
    FFTベースの高速S変換
    x: 信号（1D numpy array）
    fs: サンプリング周波数
    戻り値: (S, freqs, t)
      S: shape (n_freq, n_time) の複素数配列
      freqs: 周波数軸
      t: 時間軸
    """
    N = len(x)
    t = np.arange(N) / fs
    f_nyquist = fs / 2
    freqs = np.fft.fftfreq(N, d=1/fs)

    # FFT（周波数領域信号）
    X = np.fft.fft(x)
    
    # 出力行列
    S = np.zeros((N, N), dtype=np.complex128)

    # 各周波数でガウス窓を掛けて逆FFT
    for k, f in enumerate(freqs):
        if f == 0:
            S[k, :] = np.mean(x)  # DC成分
            continue
        # 周波数領域ガウス窓
        sigma_f = abs(f)
        g = np.exp(- (freqs - f)**2 / (2 * sigma_f**2))
        # 窓を掛けた信号を逆FFT
        S[k, :] = np.fft.ifft(X * g)
    
    # 正の周波数のみ返す
    pos_mask = freqs >= 0
    return S[pos_mask, :], freqs[pos_mask], t


# ======= デモ信号 =======
fs = 500.0
T = 2.0
N = int(T * fs)
t = np.linspace(0, T, N, endpoint=False)
# 線形チャープ + 短いバースト
f0, f1 = 20, 150
chirp = np.sin(2*np.pi*(f0*t + 0.5*(f1-f0)/T*t**2))
burst = np.zeros_like(t)
mask = (t >= 0.8) & (t < 1.0)
burst[mask] = 2.0 * np.sin(2*np.pi*80*t[mask])
x = chirp + burst

# ======= S変換 =======
S, freqs, t = fast_s_transform(x, fs)

# ======= 結果可視化 =======
plt.figure(figsize=(10, 3))
plt.plot(t, x)
plt.title("Input signal (chirp + burst)")
plt.xlabel("Time [s]")
plt.tight_layout()
plt.show()

S_mag = np.abs(S)
S_db = 20*np.log10(S_mag + 1e-10)

plt.figure(figsize=(10, 5))
extent = [t[0], t[-1], freqs[0], freqs[-1]]
plt.imshow(S_db, aspect='auto', origin='lower', extent=extent)
plt.colorbar(label='Magnitude [dB]')
plt.xlabel('Time [s]')
plt.ylabel('Frequency [Hz]')
plt.title('Fast S-Transform magnitude')
plt.tight_layout()
plt.show()

Q: S変換は、N個の時系列データを、N x N に拡張し、画像データを作成してるとも 考えられます。画像なら特徴の抽出もAIが得意とする所。心電図や脳波を解析して 異常を発見するのに利用されているんですね。

これをやってみたいのですが、生憎な事に心電図のデータなんて持ち合わして いません。でも、朝晩にサンプリングした、血圧データなら、10年分ぐらい所持 してます。余り長大だと大変なので、今年の1月から10月までのデータをS変換 したいです。時系列データを眺めると、1月ぐらいの大まかな変動が有るように 思えます。Rで時系列分析した事が有ったのですが、明確にはなりませんでした。

下記のテキストデータを取り込んで、S変換の図を表示する様に、先のサンプルを 改造してください。

sakae@lu:t$ wc blood.txt
 608  608 2432 blood.txt
sakae@lu:t$ head -n 3 blood.txt
115
106
125

データの平均値は120、標準偏差は7.9ぐらいです。特徴を余す所なく表現 できるように、元データを加工(DC成分を小さくする為、平均値を減算する等) してもかまいません。

A: 何をしたか

1. blood.txt を同じディレクトリから読み込もうとしました（見つからなかったので合成データでデモ実行）。
2. 「朝・晩サンプリング（2サンプル/日）で 2025-01-01 から 2025-10-31 まで = 608 サンプル」という想定でタイムスタンプを作成しました（あなたが示した wc の行数と合致するため、この仮定が自然と判断しました）。
3. 前処理として 平均値を引く（DC除去）＋線形トレンド除去 を行いました（ご希望どおり加工可）。
4. FFTベースの高速S変換（簡易実装）を適用し、周期（days）対時間のスペクトログラムを描画しました。
5. 図では特に「30日周期（およそ月次変動）」が見えるかどうかが分かるよう、30日線を描いています。

# s_transform_bp.py  (保存して実行)
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime, timedelta
import os

def fast_s_transform(x, fs):
    N = len(x)
    freqs = np.fft.fftfreq(N, d=1.0/fs)
    X = np.fft.fft(x)
    S = np.zeros((N, N), dtype=np.complex128)
    for k, f in enumerate(freqs):
        if abs(f) < 1e-12:
            S[k, :] = np.mean(x)
            continue
        sigma_f = abs(f)
        if sigma_f < 1e-12:
            g = np.ones_like(freqs)
        else:
            g = np.exp(- (freqs - f)**2 / (2 * sigma_f**2))
        S[k, :] = np.fft.ifft(X * g)
    pos = freqs >= 0
    return S[pos, :], freqs[pos]

# --- Main ---
fname = "blood.txt"
if os.path.exists(fname):
    data = np.loadtxt(fname).ravel()
    print(f"Loaded '{fname}' with {len(data)} samples.")
else:
    # demo data (same shape as you reported)
    N = 608
    rng = np.random.default_rng(0)
    mean_val = 120.0
    std_val = 7.9
    noise = rng.normal(0, std_val, size=N)
    samples_per_day = 2
    t_samples = np.arange(N) / samples_per_day
    monthly = 5.0 * np.sin(2 * np.pi * t_samples / 30.0)
    weekly = 2.0 * np.sin(2 * np.pi * t_samples / 7.0)
    data = mean_val + monthly + weekly + noise
    print("blood.txt not found — using synthetic demo data (mean~120, std~7.9).")

# adjust these if your timing differs
start_date = datetime(2025, 1, 1)
samples_per_day = 2

N = len(data)
timestamps = [start_date + timedelta(days=i//samples_per_day, hours=(0 if i%2==0 else 12)) for i in range(N)]

# preprocessing
x = data.astype(float)
x_mean = np.mean(x)
x_centered = x - x_mean
coefs = np.polyfit(np.arange(N), x_centered, 1)
trend = np.polyval(coefs, np.arange(N))
x_detrend = x_centered - trend

# S-transform
fs = samples_per_day
S, freqs = fast_s_transform(x_detrend, fs=fs)
periods = np.zeros_like(freqs)
with np.errstate(divide='ignore'):
    periods[1:] = 1.0 / freqs[1:]
periods[0] = np.inf

S_mag = np.abs(S)
S_db = 20 * np.log10(S_mag + 1e-12)

# plots
plt.figure(figsize=(12,3))
plt.plot(timestamps, data, label='Original')
plt.plot(timestamps, x_detrend + x_mean, label='Detrended (visual)', alpha=0.8)
plt.legend()
plt.title('Blood pressure samples (Jan–Oct 2025 assumed)')
plt.xlabel('Date')
plt.ylabel('BP')
plt.tight_layout()
plt.show()

plt.figure(figsize=(12,5))
S_db_plot = S_db[1:, :]
periods_plot = periods[1:]
im = plt.imshow(S_db_plot, aspect='auto', origin='lower',
                extent=[0, N, periods_plot[-1], periods_plot[0]])
plt.colorbar(im, label='Magnitude (dB)')
plt.ylabel('Period (days)')
plt.xlabel('Sample index (0 .. N-1)')
plt.title('S-transform magnitude (period vs time)')
if 30 >= periods_plot.min() and 30 <= periods_plot.max():
    plt.axhline(30, color='white', linestyle='--', linewidth=1)
    plt.text(5, 30*1.05, '30-day', color='white')
tick_locs = np.linspace(0, N-1, 9).astype(int)
tick_labels = [timestamps[i].strftime('%Y-%m-%d') for i in tick_locs]
plt.xticks(tick_locs, tick_labels, rotation=45)
plt.tight_layout()
plt.show()

実データが無くても雰囲気が掴める様にダミーデータを作成してくれている。 親切なこったい。ダミーの作り方が面白いなあ。月の変動、週の変動と多少のノイズ を加算してる。毎日が日曜日の年金生活者は、月締めのストレスにも週報の面倒 くささも無縁なんだけどな。サラリーマンとかだと、こういう社会性ストレスにも 影響しますよって仮定なんだな。

下記の様に呼出を変更すると、ウェーブレット版になる。

S, freqs = wavelet_transform(x_detrend, fs=fs, n_freqs=200, omega0=6.0)

def wavelet_transform(x, fs, n_freqs=200, f_min=None, f_max=None, omega0=6.0):
    """
    FFTベースの連続ウェーブレット風変換（Morlet-like）。
    - x: 1D real signal (length N)
    - fs: sampling rate (samples per unit time, e.g. samples/day)
    - n_freqs: 出力する周波数ビン数（正の周波数のみ）
    - f_min, f_max: 出力周波数レンジ（cycles per unit time）。指定しなければ [1/N*fs, fs/2]。
    - omega0: モルレ母関数の中心角周波数（典型値 6.0）。大きいほど周波数分解能↑、時間分解能↓。
    戻り値:
    - S: shape (n_freqs, N) complex array (wavelet coefficients)
    - freqs: length-n_freqs array of positive frequencies (cycles per unit time)
    """
    x = np.asarray(x, dtype=float)
    N = x.size
    # FFT of signal
    X = np.fft.fft(x)

    # default frequency range
    nyq = fs / 2.0
    df_min = fs / N
    if f_min is None:
        f_min = max(df_min, 1e-12)
    if f_max is None:
        f_max = nyq

    # frequency vector (positive frequencies)
    freqs = np.linspace(f_min, f_max, n_freqs)

    # time vector centered at zero for wavelet construction (length N)
    t = (np.arange(N) - N//2) / fs  # seconds (or days) centered at 0

    # precompute FFT frequencies for wavelet FFTs (same length N)
    # We'll compute FFT(psi) by fftshifting the time-domain psi constructed on t centered at zero
    S = np.zeros((n_freqs, N), dtype=np.complex128)

    # For a Morlet-like wavelet: psi(t) = C * exp(i*omega0*t/s) * exp(-t^2/(2*(s*sigma0)^2))
    # We choose scale s so that center frequency f = omega0 / (2*pi*s)  => s = omega0/(2*pi*f)
    # Then time-width sigma_t = s (here we use s as time-scale modulated by a factor)
    # Practical implementation: for target frequency f, set scale s = omega0/(2*pi*f)
    # and sigma = s (or s/np.sqrt(2)); we choose sigma = s/ np.sqrt(2) to have reasonable bandwidth.
    for i, f in enumerate(freqs):
        if f <= 0:
            S[i, :] = 0.0
            continue
        # scale corresponding to target frequency
        s = omega0 / (2.0 * np.pi * f)
        sigma = s / np.sqrt(2.0)  # Gaussian std in time
        # construct wavelet centered at t=0 (length N)
        # normalization: 1/sqrt(s) typical for CWT energy normalization
        norm = 1.0 / np.sqrt(s)
        psi_t = norm * np.exp(-0.5 * (t / sigma)**2) * np.exp(1j * 2.0 * np.pi * f * t)
        # shift psi so its time-zero is at index 0 for convolution via FFT:
        psi_shifted = np.fft.ifftshift(psi_t)  # move center to array index 0.. for proper FFT conv
        Psi = np.fft.fft(psi_shifted)
        # convolution x * conj(psi)  => in frequency domain: X * conj(Psi)
        coeff = np.fft.ifft(X * np.conjugate(Psi))
        S[i, :] = coeff

    return S, freqs

ウェーブレット変換は、無線機のバンドモニターに採用されてるな。貧乏な 無線家は、せいぜいSDRのウォーターフォールを眺めるぐらいが関の山だ。 もっと一般的には、チャッピー君の音声認識に人知れず使われている。

Q: S変換とは直接関係無い事ですが、質問します。

サンプリング定理によれば、原信号の成分の最大周波数 fmax の2倍（2fmax）よりも 高い周波数で標本化した信号は、元の信号に復元できるとあります。 当然この仮定として、サンプリング周期は等間隔だと理解してます。

今回の血圧データは、2回/日でサンプリングしてます。実際の採取時刻は、朝4時と 夕21時です。日に対して、サンプリング周期が等間隔ではありません。 この様なデータから、日のデータを正しく復元できるものなんでしょうか？

こういう事例は、工業界では無いと思われますが、あえて問う次第です。

A: とても良い質問です。 おっしゃる通り、今回のような 非等間隔サンプリング（4時と21時） の場合、 「サンプリング定理（ナイキスト＝シャノンの標本化定理）」は /そのままでは成立しません/。 ただし、「復元がまったく不可能」というわけでもなく、いくつかの観点があります。

---

## 🧩 1. サンプリング定理の前提とその破れ方

標本化定理が保証しているのは次の条件です：

周期的に等間隔 Δt でサンプリングされた信号 x[n] があるとき、 もし元信号の周波数成分が fmax < 1/(2Δt) に完全に制限されているなら、 x[n] から元の信号を完全に再構成できる。

つまり：

• サンプリング周期は /一定（等間隔）/。
• 信号は /バンドリミット（高周波成分が有限）/。

このうち 非等間隔サンプリング は明確に前提を破っています。 この場合、通常のFFTやナイキスト理論による完全再構成は保証されません。

---

## 📉 2. 今回の例：1日2回（4時・21時）

### サンプリング時刻

• 1点目：朝4時 → 時刻=4/24=0.1667日
• 2点目：夜21時 → 時刻=21/24=0.875日
• 次の日の朝：+1日で 1.1667日
• したがって周期が「0.7083日」と「0.2917日」で交互になります。

つまり、/ほぼ2サンプル/日/ ですが、/周期が交互に違う/。

これを等間隔 0.5日のサンプリングに近似するとき、 最大で ±0.2083日のタイミング誤差（約5時間）が入ることになります。 このズレが「擬似的なジッター」として作用し、 スペクトル解析（特に高周波側）に/歪み（スペクトル拡散・偽ピーク）/を生じさせます。

---

## 🧮 3. 非等間隔サンプリング信号の取り扱い方法

この問題は 天文観測データ/・/生理時系列（脈拍・血糖）/・/金融データ など、 工業以外の分野では頻繁に出ます。 次のようなアプローチがあります。

### (a) 再サンプリング（等間隔補間）

• 時刻データ（例：=t = [0.1667, 0.875, 1.1667, 1.875, …]=）を明示し、 スプライン補間などで等間隔格子上に補間する。
• 得られた等間隔信号に対してFFT・S変換・ウェーブレット変換などを適用。
• 短期的な周期解析（例：7日周期、30日周期）にはこの方法で十分有効。

⚠️ 注意：補間により高周波成分が人工的に滑らかになるため、 /時間分解能は実際よりも劣化/します。

---

### (b) 非等間隔フーリエ変換（Nonuniform FFT; NFFT）

• サンプリングが非等間隔でも、フーリエ成分を直接求めるアルゴリズムがあります。
• NumPyには標準で入っていませんが、 pynfft や nfft が利用可能です。
• これを使うと補間なしで周波数解析可能ですが、 実装と理論が少し重い（S変換・ウェーブレットに組み込むには拡張が必要）。

---

### (c) 非等間隔小波変換

• Morletウェーブレットなどは、時刻 t を明示的に使って構築できるため、 サンプリング点が不均等でも離散積分の形で計算可能です。
• ただし、NumPyで簡単に実装するのは手間（=Δt=が変化するたびに重み補正が必要）。

---

## 📊 4. 現実的な結論（今回の血圧データの場合）

目的	対応策	備考
----------------	----------------------------	----------
周期性（例：月次・週次）を見たい	スプライン補間→等間隔化→FFT/S変換/ウェーブレット	✅ 現実的・十分正確
ミリ秒単位の高周波解析	不可（サンプリング定理以前にサンプリング不足）	—
正確な位相再構成・原波形再現	不可能（サンプリング定理の破れ）	—
傾向・周期推定・異常検知	問題なし（周期 > 2日程度なら誤差小）	✅

---

## 🔍 5. まとめ（要点）

• サンプリング定理は 等間隔サンプリング が前提。
• 4時と21時という非等間隔は「軽度のサンプリングジッター」と見なせる。
• 週単位・月単位など/長周期の変動解析には実用上問題ない/。
• FFT/S変換を使うときは、一度 補間して等間隔化 してから解析するのがベター。

今回取り上げたS変換は、やってる事が良く解るよ。入力されたデータを一度FFTして 周波数ドメインなデータに変換する。それをそのまま逆FFTすれば元に戻る。

これじゃ意味が無いので、FFTした周波数ドメインなデータに窓関数を乗算する。 窓関数は、周波数に対応する特性を持つフィルターだ。で、この結果を逆FFT。 これを周波数の隅々まで適用する。

このS変換の論文を取り上げた博士が、悔しがっていたぞ。誰でも考え付きそうだけど 今迄無かった方が不思議ってね。まあ、自分用の道具箱の道具を一つ増やした。 有益だったんで、論文にして公開するよってOSSの思想なんだな。

そうか、有益なライブラリィーとかは、ただGithubに掲載するだけじゃなくて、 論文にすればいいのか。論文が採択されれば箔が付くもんな。数多ある食堂でも ミシュランの星を獲得すれば、関係者のお墨付きって訳だな。博士は何かと忙しい から、何か作って提供したら、良い事があるかも。

オイラーも昔、同様な 思想で、ローパスフィルターを作った事が有る。 一度FFTする。その結果で不要な高域周波数成分にZEROを書き込んでしまう。そうして 逆変換すれば、フィルターされた事になる。論文にしたら採択されたかも。