ラズパイのGPIO割り込み検出とsystemd設定について

ラズパイにseismicサービスを組み込む時に関連したメモ

① 現状Rust(rppal)ではGPIOで割り込みを検出する手段は提供されていない様子

作ればいいんだろうけど、今の所クレートは存在していないから、従来通りそこだけはPythonのサービスを起動、以下のソースで個別にサービス定義して起動時に実行させておく

#
# wait switch push interrupt and issue shutdown command
#
import time
import datetime
import os
import sys
import RPi.GPIO as GPIO
import subprocess

# Shut down sw is assigned to GPIO17
# GPIO initialize
SHUTDOWN = 17

GPIO.setwarnings(False)
GPIO.setmode(GPIO.BCM)
GPIO.setup(SHUTDOWN, GPIO.IN)
GPIO.setup(SHUTDOWN, GPIO.IN, pull_up_down=GPIO.PUD_UP)

def handle_sw_input():
# wait key inout event
    def switch_callback(gpio_pin):
        subprocess.call('sudo shutdown -h now', shell=True)
#
    GPIO.add_event_detect(SHUTDOWN, GPIO.FALLING,bouncetime=500)
    # when the sw was pushed, call the 'call back routine' 
    GPIO.add_event_callback(SHUTDOWN, switch_callback) 
    return

handle_sw_input()

while True:
    time.sleep(100)

② 今更ですが。/lib/systemd/systemのserviceファイルを変更した時には、以下の処理が必要(既存ファイルの書き換え替えが反映されなかった)

・サービスを登録する(編集したときの反映時にも必要)

$ sudo systemctl daemon-reload

$ sudo systemctl enable seismic.service

 

admin

表示まで含めてRustで実装

ssd1306への表示データをchannel使って別スレッドに送って、その別スレッドでssd1306に震度データを表示させるような構成にして、Python版からRust版に置き換え完了

 

ただし、シャットダウンスイッチを扱うコードはまだ作成していないので、それを追加しないとスタンドアロンで運用はできない

https://github.com/chateight/seismic/blob/main/Rust_version/main_disp.rs

Cargo.tomlの中身は、

[package]
name = "seismic_refatoring"
version = "0.1.0"
edition = "2021"

[dependencies]
chrono = "0.4.38"
rppal = "0.19.0"
linux-embedded-hal = "0.4.0"
embedded-graphics = "0.8.1"
machine-ip = "0.2.1"
ssd1306 = "0.9.0"

震度計プロジェクトも概ね終わりかな

 

admin

 

今の簡易震度計算ロジックの精度の検証

Rustで記述したPython版の震度計算ロジックの精度検証のために、公開データを使って計算させてみた

https://www.data.jma.go.jp/eqev/data/kyoshin/jishin/2401011610_noto/index.html

ここから波形データをcsv形式でダウンロードできるので、外部ファイル読み込みでデータを取り込ませて計算させた

コードは以下のリンクから、

https://github.com/chateight/seismic/tree/main/Rust_version

気象庁の参考震度情報と比較してもほぼ一緒と言えるレベルなので、実用上は問題ないと言えます

あえて修正が必要と思えるのは、

① 加速度センサーと震度波形データのスケール合わせ

② 実際のハードでは加速度計の出力におよそ30Hzの一次ローパスフィルタが入っているので実際の波形データを前処理してみるか

ぐらいですが、まあみたところそこまで必要のなさそうなレベルです、実際の地震計は設置だけでも大変なレベルなのに、それを自宅の部屋で使おうというぐらいだから

P.S. 2024/11/14

微妙に違和感あったので、見直すとVectorでFIFOで使っているのにわざわざリングバッファ扱いでポインタ使って正しくない場所からadc_data持って来ていたので修正

 

admin

震度計のPythonコードをRustに書き換え

ラズパイzeroで作った震度計のコードをRustで書き換えてみた、

https://github.com/p2pquake/rpi-seismometer/blob/master/seismic_scale.py

とりあえず変更点は何もなく置き換えたつもり、

https://github.com/chateight/seismic/tree/main/Rust_version

x/y/zに200.0の固定値を与え続けて、M1 Macで動作させた実行初期の出力される値はこんな感じで、小数点以下二桁目から値が違うのはまだロジックが完全に一致していない可能性がある、実用的には何の問題もないのだけれども

% python -u "/Users/usamiryuuichi/python/seismic.py"
2024-11-08 19:49:45.068661 scale: 0  frame: 0
2024-11-08 19:49:45.169671 scale: 0  frame: 20
2024-11-08 19:49:45.269686 scale: 0  frame: 40
2024-11-08 19:49:45.369185 scale: 4.141110582455145  frame: 60
2024-11-08 19:49:45.469702 scale: 5.676822095241448  frame: 80
2024-11-08 19:49:45.569745 scale: 5.676822095241448  frame: 100
2024-11-08 19:49:45.669175 scale: 5.676822095241448  frame: 120
2024-11-08 19:49:45.769694 scale: 5.676822095241448  frame: 140


% cargo run --release
    Finished `release` profile [optimized] target(s) in 0.01s
     Running `target/release/seismic`
time: 2024-11-08 22:28:28.151209 UTC, scale: 0, frame: 0
time: 2024-11-08 22:28:28.252217 UTC, scale: 0, frame: 20
time: 2024-11-08 22:28:28.351601 UTC, scale: 0, frame: 40
time: 2024-11-08 22:28:28.451895 UTC, scale: 4.18249, frame: 60
time: 2024-11-08 22:28:28.552350 UTC, scale: 5.681231, frame: 80
time: 2024-11-08 22:28:28.652210 UTC, scale: 5.681231, frame: 100
time: 2024-11-08 22:28:28.751913 UTC, scale: 5.681231, frame: 120
time: 2024-11-08 22:28:28.852249 UTC, scale: 5.681231, frame: 140

次のステップの能登の地震波形から、気象庁の震度データ波形と簡易計算結果とを比較してみること、簡易方式は30Hzあたりでアナログ的に一次のLPF掛けてるから、波形データに対して同等の前処理が必要かもしれないし、近似度がそれほでないデータならそれはおそらく無意味

 

admin

別のスレッドにメッセージを送る(@Rust)

現在のRustは非同期処理は得意な分野で(OS開発にも使われるぐらいだから当然)、非同期処理だけで書籍があるぐらい奥も深い

ここでは単純に別のスレッドにメッセージを送ることだけをやってみた

Rustではスレッド間通信はクロージャーを使うようで、rust-lang.orgのサンプルコードを多少変えています

//
// https://doc.rust-lang.org/std/thread/fn.spawn.html
//
use std::thread;
use std::sync::mpsc::channel;
use std::time::Duration;

fn main() {
    // チャネルを作成
    let (tx, rx) = channel();

    // レシーバースレッドをspawn
    let _receiver = thread::spawn(move || {
        let value: String = rx.recv().expect("Unable to receive from channel");
        rcv(value.clone());
    });

    // センダースレッドをspawn
    let _sender = thread::spawn(move || {
        tx.send("Hello, thread".to_owned())
            .expect("Unable to send on channel");
    });

    // センダースレッドの終了を待つ
    //sender.join().expect("The sender thread has panicked");

    // レシーバースレッドの終了を待つ(但し、rcv関数の終了を待たない)
    //receiver.join().expect("The receiver thread has panicked");

    println!("reach to the end");
    thread::sleep(Duration::from_secs_f32(1.5));
}

fn rcv(value: String) {
    thread::sleep(Duration::from_secs_f32(1.0));
    println!("recieved : {}", value);
}

このコードではsenderもreceiverも終了を待たないようにしていますが、その場合にはrcv関数の方が時間が遥かにかかるので、main関数の最後の時間待ちがないとrev()の実行は途中で打ち切られます

今回のケースでは相互にデータをやり取りや完了を待つわけではなく、main関数からrcv関数にメッセージを送りたいだけなのでこのような構成で十分です

さらに言えば、senderスレッドを別に起こさなくとも直接tx.send()を実行しても構いません

もし完全に非同期処理が必要ならばasync/awaitの出番ですが、これはgolangと結構類似しているように思います

 

admin

Rustでspi経由でmcp3204からのデータを読み取る

最初このクレートを使おうと思いましたが、

https://github.com/trashware/mcp3xxx-rs

開発が5年前で、かなり古くて依存先のrppalの版数が0.11.0までになってます、rppalではspi経由のデータやり取りはコマンドストリームの送信と結果のストリームが全二重で処理できるから、別に直接ビット列を用意すればクレート使わなくても良いわけで、最終的に以下のコードで読み取りできました

Cargo.tomlではrppalは最新版を指定しています

use rppal::spi::{Bus, Mode, SlaveSelect, Spi};
use std::thread;
use std::time::Duration;

const CHANNEL_0: u8 = 0x06 | 0x00; // Command byte for single-ended mode

fn main() {
    // Create a SPI object
    let spi = Spi::new(Bus::Spi0, SlaveSelect::Ss0, 1_000_000, Mode::Mode0).unwrap();

    // Define the buffers for reading and writing
    let mut read_buffer = [0u8; 3];
    let write_buffer = [CHANNEL_0, 0x00, 0x00];

    // Poll the sensor every n seconds
    let poll_rate = Duration::from_secs(1);
    loop {
        // Perform the SPI transfer
        spi.transfer(&mut read_buffer, &write_buffer).unwrap();

        // Extract the 12-bit ADC value
        // Note: MCP3204 returns 12-bit data, but the second byte contains status bit(0) and the MSB of the result.
        let msb = (read_buffer[1] & 0x0F) as u16;
        let lsb = read_buffer[2] as u16; // Convert third byte to u16
        let result = (msb << 8) | lsb;

        println!("Sensor value: {}", result);      

        println!("{:08b}, {:08b}", read_buffer[1], read_buffer[2]);

        thread::sleep(poll_rate);
    }
}
[dependencies]
rppal = "0.19.0"

spi上のストリーム情報はデバイスの仕様から、SGLモードを使うので1バイト目にはstartビットと合わせて6を指定すれば良いことになります

コードサンプルではch0の読み取りを行なってますが、このビット列を変更することで他のチャネル情報も取得できます

以下は実行した時の様子

$ ./adc_rppal
Sensor value: 2867
00001011, 00110011
Sensor value: 2867
00001011, 00110011
Sensor value: 2867
00001011, 00110011
Sensor value: 2867
00001011, 00110011
Sensor value: 2868
00001011, 00110100

 

admin

VMwareのUbuntuがクラッシュした

Intel MacのVMwareのUbuntuが立ち上がらなくなった、「コピーしますか」とかいうポップアップ出ていてそれが関連するかどうかだけれども、ともかくも身動き取れなくなったからSSDにバックアップしていた仮想マシンバックアップをコピー戻してで修復、区画は大きめに取っていても実際のサイズでコーピーが終了するようだ、つまり全領域の復元はやらないということ

ついでにVMwareのアップデートも適用、Ubuntu20系のサポートは来年5月までだからそれほど残り時間はない

 

admin

能登地震の波形のスペクトルを求めてみた

正弦波の合成だけでは実用性はないので、今年の元旦の能登地震の公開されている波形データからスペクトルを求めてみた

公開先はこちらですが、トップにある輪島市のデータを使っています、形式はcsvなのでヘッダー情報除けばPythonやRustで簡単に処理できます

・ヘッダー情報

SITE CODE= 67016輪島市門前町走出        ,37.4962,137.2705,15.86,7.6,45292.67387,8.93,15.21
 LAT.=   37.2871,,,,,,,
 LON.=  136.7680,,,,,,,
 SAMPLING RATE= 100Hz,,,,,,,
 UNIT  = gal(cm/s/s),,,,,,,
INITIAL TIME = 2024 01 01 16 10 10,,,,,,,
 NS,EW,UD,,,,,

 

気象庁|強震観測データ|2024/1/1 石川県能登地方の地震

Pythonで全時間領域を対象にしてみたもの(画像はUD軸)

リンク先にある波形と比較するとほぼ類似の形状になってます

 

以下はRustで部分切り出し(全体で100サンプル/secで120,000フレーム、つまり120秒分のデータがありますが一番のピークの3,000~4012部分の切り出し)を、窓処理してFFTしてみたもの

<NS波形>

<NSスペクトラム>

以下はEW

以下はUD

 

NS/EWに比較すると周波数成分が高い方に出てきます、Pythonの全時間軸ともちろん傾向的には同じようなスペクトルになっています

X軸の20がほぼ10Hzに相当します

 

以下はPythonのコードになりますが、対話形式でほぼGeminiで作成させたコードがそのまま動きました、きちんと境界条件を与えてやれば人間がコードを書く手間を大幅に省力化できます、コードは読めないとダメですが

import numpy as np
import matplotlib.pyplot as plt

def fft_csv(filename, column_index, sampling_rate):
  """
  CSVファイルの特定列データをFFTする関数

  Args:
    filename: CSVファイルのパス
    column_index: FFTしたい列のインデックス (0から始まる)
    sampling_rate: サンプリングレート

  Returns:
    freqs: 周波数
    fft_result: FFTの結果
  """

  # CSVファイルを読み込む
  data = np.loadtxt(filename, delimiter=',', skiprows=1)  # ヘッダー行をスキップ

  # 指定した列のデータを取得
  target_data = data[:, column_index]

  # FFTを実行
  fft_result = np.fft.fft(target_data)

  # 周波数軸を作成
  N = len(target_data)
  freqs = np.fft.fftfreq(N, 1.0/sampling_rate)

  return freqs, fft_result

# 使用例

if __name__ == "__main__":
    filename = "noto.csv"  # CSVファイル名
    column_index = 2  # FFTしたい列 (3列目) 0 :NS/1 :EW/2 :UD
    sampling_rate = 100  # サンプリングレート

    freqs, fft_result = fft_csv(filename, column_index, sampling_rate)

# 結果をプロット (両対数プロット)
    plt.loglog(freqs, np.abs(fft_result))
    plt.xlabel("Frequency [Hz]")
    plt.ylabel("Amplitude")
    plt.title("FFT of Column {}".format(column_index+1))
    plt.grid(True)
    plt.show()

 

admin

 

スペクトルグラフをLogスケールにしてみる

普通は周波数スペクトルのグラフはLogスケールなのでそれでみてみる

ソースコードの変更箇所は、前回のコードから以下の一箇所だけ変更で、10を底にする普通のLogスケールに変換しています

    // absolute value calc
    let y: Vec<f32> = buffer.iter().map(|z| z.norm().log(10.0)*20.0).collect();

<結果のグラフ>

信号レベルからかなり低いのですが、そこそこのレベルで『ゼロ』レベルが推移していますが、これはf32であることによる丸め誤差で、f64にすると景色が変わります、実用上は検出した信号に対して100db以上のダイナミックマージンあればこれ以上の精度は不要だろうし、計算速度の点からもデメリットがあると思う

 

 

admin

rustfftで窓をかける

実際のFFTでは無限に継続する波形というのはないので、実用的なアプリケーションでは信号の切り出しが必要です

ところでrustfftには窓機能のメソッドはないので、個別に実装が必要ですが、一番ポピュラーと思うhanning窓を適用してみました

hanning窓は(1-cos(x))/2で表現されますが、前回の正弦波の合成波形にこの窓を適用します

<窓適用後の入力波形>

<FFTの結果>

比較用の窓処理なしのスペクトル

窓をかけるとピークレベルはほぼ半分になります(対称波形のエンベロープの面積から)、また多少スペクトルの根元が広がっています、予測ではもっと広がると思ってた

P.S. 2024/10/28

スペクトルの分散の考え方は、このケースでは正弦波を窓関数でAM変調していると捉えれば、キャリア周波数の両サイドに山がいくつか見えるはずで、周波数分解能が足りないからダンゴになって見えているだけです

<コード>

正弦波の合成のところで窓処理も同時に行ってます

use plotters::prelude::*;
use rustfft::{num_complex::Complex, FftPlanner};
use std::f32::consts::PI;
use std::time::Instant;

const N: usize = 256; // FFT size
const SAMPLE_RATE: f32 = 256.0; // sampling rate
const INPUT_FREQ: f32 = 20.0; // input frequency

fn complex_vector_with_hanning(length: usize, dt: f32, frequency: f32) -> Vec<Complex<f32>> {
    (0..length)
        .map(|i| {
            let t = i as f32 * dt;
            let phase = frequency * 2.0 * PI * t;
            let phase2 = 3.0 * frequency * 2.0 * PI * t;
            Complex::new(
                (f32::sin(phase) + f32::sin(phase2)) * 0.5 * (1.0 - f32::cos(2.0 * PI * t)),
                0.0,
            )
        })
        .collect()
}

fn re_vector(complex_numbers: Vec<Complex<f32>>) -> Vec<f32> {
    complex_numbers.iter().map(|z| z.re).collect()
}

fn draw(
    x: Vec<usize>,
    y: Vec<f32>,
    f_name: &str,
    cap: &str,
) -> Result<(), Box<dyn std::error::Error>> {
    let image_width = 1080;
    let image_height = 720;

    let root = BitMapBackend::new(f_name, (image_width, image_height)).into_drawing_area();

    root.fill(&WHITE)?;

    //   https://qiita.com/lo48576/items/343ca40a03c3b86b67cb
    let (y_min, y_max) = y
        .iter()
        .fold((0.0 / 0.0, 0.0 / 0.0), |(m, n), v| (v.min(m), v.max(n)));

    let caption = cap;
    let font = (“sans-serif”, 20);

    let mut chart = ChartBuilder::on(&root)
        .caption(caption, font.into_font())
        .margin(10)
        .x_label_area_size(16)
        .y_label_area_size(42)
        .build_cartesian_2d(*x.first().unwrap()..*x.last().unwrap(), y_min..y_max)?;

    chart.configure_mesh().draw()?;

    let line_series = LineSeries::new(x.iter().zip(y.iter()).map(|(x, y)| (*x, *y)), &RED);
    chart.draw_series(line_series)?;

    Ok(())
}

fn main() {
    // fft data preparation with window
    let sr = 1.0 / SAMPLE_RATE;

    let mut buffer = complex_vector_with_hanning(N, sr, INPUT_FREQ);
 
    // draw input signal
    let buf_re = re_vector(buffer.clone());

    let n = N;
    let f_input = “plot_input.png”;
    let cap = “fft input”;
    let x: Vec<usize> = (1..=n).collect();

    let _ = draw(x, buf_re, &f_input, &cap);

    // fft mode setting
    let mut planner: FftPlanner<f32> = FftPlanner::new();
    let fft = planner.plan_fft_forward(N);

    let start_time = Instant::now();
    
    // fft exec
    fft.process(&mut buffer);

    let end_time = Instant::now();
    println!(“Elapsed time for FFT: {:?}”, end_time.duration_since(start_time));

    // absolute value calc
    let y: Vec<f32> = buffer.iter().map(|z| z.norm()).collect();

    // draw fft graph output
    let n = N;
    let file_fft = “plot_fft_hanning.png”;
    let cap = “fft hanning result”;
    let x: Vec<usize> = (1..=n).collect();
    let _ = draw(x, y, &file_fft, &cap);

}

 

admin