FFTフィルタではなく、IIRフィルタを使った場合のLPFを試してみた。
次数が大きくなるほど位相遅れが大きくなり、下記では12次までやってみたが、14次は計算が破綻した（位相が１周ひっくり返るからかな）
例えば２次だと普通に２次のLPFですね（下はFFT後の縦軸がdBになってないです）
８次のフィルタだと位相は反転しているくらいの遅れが見られますね。ふむふむ。
FFTして周波数帯域をカットしてからの逆FFTで信号を戻す作戦。
単純な波形においては問題がなさそうである。
3倍高調波である75Hzまで入れると・・・（カットは100Hzで）
25Hzも75Hzも大丈夫そうな感じではある。
pythonは配列の要素にアクセスせず、配列のまま計算できるのが良い所である。
import numpy as np
import matplotlib.pyplot as plt
freq = 1.0 #周波数として
time=np.arange(0.000,0.001*2**10,0.001) #時間の配列を作る
Signal = np.sin(2 * np.pi * freq * time) #正弦波を作る
plt.plot(time,Signal)
plt.show()
次のような配列要素ごとの計算をしなくて良いので、記述がやや楽。
for i = 0 to Ubound(time())
Signal(i) = sin(2 * pi() * freq * time(i))
next i
正弦波は単一周波数しか含まれないので、矩形波や三角波を作ろうと思うと・・・・
import numpy as np
import matplotlib.pyplot as plt
freq = 2.0 #周波数として
time=np.arange(0.000,0.001*2**10,0.001) #時間の配列を作る
Signal = (time * freq) % 1 #余りを求める
plt.plot(time,Signal)
plt.show()
これを0.5ずらして絶対値とったら三角波ですかね。
import numpy as np
import matplotlib.pyplot as plt
freq = 2.0 #周波数として
time=np.arange(0.000,0.001*2**11,0.001) #時間の配列を作る
Signal = 1-np.abs((time * freq/2) % 1-0.5)*2
plt.plot(time,Signal)
plt.show()
のこぎり波など不連続データが多いと周波数帯域が非常に高い所まで必要になるので、一旦三角波を使ってフィルタ特性を見てみようか、な。
検証はFFTをするためにデータサイズを2のべき乗の数にしているが、実際の信号はそんなに都合よくデータ数がないので、最後はFIRフィルタとかIIRフィルタとかを作って連続信号に対して処理できればいいかな、なんて思い始めている。
入力信号は25Hzで振幅1の正弦波に60Hzで振幅0.5の正弦波を加算したものである。
SignalF = 25.0
SignalA = 1.0
NoiseF = 60.0
NoiseA = 0.5
time=np.arange(0.000,0.001*2**12,0.001)
SignalSUM = SignalA * np.sin(2*np.pi*SignalF*time) + NoiseA * np.sin(2*np.pi*NoiseF*time)
左上の入力信号をFFTして、左下のようなグラフを得る。
左下のグラフの60Hzの山をカットして右下の周波数特性を得る。
その信号を逆FFTして、波形に戻すとローパスされた信号が得られる。
ノイズの信号を大きくしていったらどうなるのかなと思って信号振幅を変えてみた。
信号：ノイズ＝１：１
信号：ノイズ＝１：１０、なんとかなる。
信号：ノイズ＝１：１００ FFTで一律ゼロにする部分がうまく行かないので信号が歪んでくる。
40Hzでスパッと切っているのだが、40Hzの階段が40Hzのノイズを与えてしまうので、しょうがない。
import numpy as np
import matplotlib.pyplot as plt
#from scipy import signal
import csv
SignalF = 25.0
SignalA = 1.0
NoiseF = 60.0
NoiseA = 100.0
GAmp = SignalA+ NoiseA
time=np.arange(0.000,0.001*2**12,0.001)
SignalSUM = SignalA * np.sin(2*np.pi*SignalF*time) + NoiseA * np.sin(2*np.pi*NoiseF*time)
#plt.plot(time,SignalSUM, '-')
F=np.fft.fft(SignalSUM)
freq = np.linspace(0, 1000, len(time))
Fc= int(40 * len(time)/1000)
Fmod = F.copy()
Fmod[Fc:(len(time)-Fc)]=0+0j
#plt.plot(freq,F)
#plt.xlim(0,80)
G=np.fft.ifft(Fmod)
ax1 = plt.subplot(2,2,1)
ax1.plot(time, SignalSUM, "-")
ax2 = plt.subplot(2,2,2)
ax2.plot(time, G , "-")
ax1.set_xlim(0,0.5)
ax2.set_xlim(0,0.5)
ax1.set_ylim(-GAmp,GAmp)
ax2.set_ylim(-GAmp/10,GAmp/10)
ax3=plt.subplot(2,2,3)
ax3.plot(freq, F)
ax3.set_xlim(0,100)
#ax3.set_ylim(-100,100)
ax4=plt.subplot(2,2,4)
ax4.plot(freq, Fmod)
ax4.set_xlim(0,100)
#ax4.set_ylim(-10,10)
print(len(time))
plt.show()
MathJaxというものがある。TeXで数式をキレイに書いてくれるものである。
MathJaxを読み込んで終わりではあまり面白くないので、tex形式で入力した式を書いてくれるグッズを作ってみた。ChatGPTで。
仕事でちょっとほしかったのだが、悩んでいてそうだ、こういうのこそChatGPTだとおもって作ってもらったらあっさりできた。
〜〜〜以下HTMLコード〜〜〜〜
<!DOCTYPE html>
<html>
<head>
<title>MathJax TeX Test Page</title>
<script src="https://polyfill.io/v3/polyfill.min.js?features=es6"></script>
<script id="MathJax-script" async src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>
</head>
<body>
<textarea id="math" rows="4" cols="50"></textarea>
<button onclick="renderMath()">Render</button>
<script>
function renderMath() {
var math = document.getElementById('math').value;
var output = document.getElementById('output');
output.innerHTML = '';
MathJax.texReset();
var options = MathJax.getMetricsFor(output);
options.display = true;
MathJax.tex2chtmlPromise(math, options).then(function (node) {
output.appendChild(node);
MathJax.startup.document.clear();
MathJax.startup.document.updateDocument();
});
}
</script>
</body>
</html>
シリアル信号からUSB HIDデバイスに信号を変換するICというものがすでに存在していた。
秋月電子に。
モジュール状態ではAliexpressからも入手可能なようである。
しかし、9600baudがデフォルトなので、今使いたい天秤だと設定変更が必要なので、マイコンで作っていくのが良さそうかな。なんて。
Savitzky-Golayフィルタ関数を定義します。この関数では、入力データの配列、窓幅、多項式次数、微分の次数を引数として取ります。
Function SavitzkyGolayFilter(data As Variant, window As Integer, order As Integer, deriv As Integer) As Variant
Dim i As Integer, j As Integer, k As Integer
Dim sum As Double, a As Double, b As Double
Dim coef() As Double
Dim filteredData() As Double
Dim n As Integer
n = UBound(data)
ReDim filteredData(n)
ReDim coef(order + 1, window)
For i = 0 To n
sum = 0
For j = -window To window
If (i + j >= 0 And i + j <= n) Then
For k = 0 To order
coef(k, j + window) = coef(k, j + window) + (j ^ k)
Next k
sum = sum + data(i + j)
End If
Next j
filteredData(i) = sum
Next i
For i = 0 To order
For j = 0 To window
coef(i, j) = coef(i, j) / coef(0, window)
Next j
Next i
For i = 0 To n
a = 0
For j = -window To window
If (i + j >= 0 And i + j <= n) Then
a = a + coef(deriv, j + window) * data(i + j)
End If
Next j
filteredData(i) = a
Next i
SavitzkyGolayFilter = filteredData
End Function
実際にフィルタリングを行うには、以下のように関数を呼び出します
Sub ApplySavitzkyGolayFilter()
Dim inputData() As Double
inputData = Range("A1:A100").Value
Dim window As Integer
window = 5
Dim order As Integer
order = 2
Dim deriv As Integer
deriv = 0
Dim outputData() As Double
outputData = SavitzkyGolayFilter(inputData, window, order, deriv)
Range("B1:B100").Value = outputData
End Sub
Excelでは、NORM.INV関数を使用して、指定された平均値と標準偏差に基づいて正規分布する乱数を生成することができます。以下の手順に従って、正規分布する乱数を作成してください。
-
平均値と標準偏差を決定します。例えば、平均値を10、標準偏差を2に設定したい場合は、それぞれセルに入力します。
-
NORM.INV関数を使用して、正規分布する乱数を生成します。関数の構文は以下の通りです。NORM.INV(probability, mean, standard_dev)probability引数は、0から1までの確率値です。mean引数は、正規分布の平均値を指定します。standard_dev引数は、正規分布の標準偏差を指定します。例えば、セルA1に
=NORM.INV(RAND(),10,2)と入力することで、セルA1に平均値10、標準偏差2の正規分布する乱数を生成することができます。 -
RAND()関数を使用して、0から1までのランダムな数値を生成します。この数値は、NORM.INV関数のprobability引数に渡されます。RAND()関数を使用することで、平均値と標準偏差に従う正規分布する乱数を生成することができます。例えば、セルA1に
=NORM.INV(RAND(),10,2)と入力することで、セルA1に平均値10、標準偏差2の正規分布する乱数を生成することができます。 -
上記の手順を繰り返すことで、必要な数の正規分布する乱数を生成することができます。例えば、10個の正規分布する乱数を生成する場合は、セルA1に上記の式を入力し、セルA2からA10に同じ式をコピー＆ペーストします。
以上の手順に従うことで、Excelで正規分布する乱数を作成することができます。
サブルーチンとして次を定義しておく。
Sub QuickSort(arr As Variant, firstIndex As Long, lastIndex As Long)
Dim pivot As Variant
Dim pivotIndex As Long
Dim leftIndex As Long
Dim rightIndex As Long
If firstIndex < lastIndex Then
pivot = arr(lastIndex)
pivotIndex = firstIndex - 1
leftIndex = firstIndex
For rightIndex = firstIndex To lastIndex - 1
If arr(rightIndex) <= pivot Then
pivotIndex = pivotIndex + 1
Swap arr(pivotIndex), arr(rightIndex)
End If
Next rightIndex
pivotIndex = pivotIndex + 1
Swap arr(pivotIndex), arr(lastIndex)
Call QuickSort(arr, firstIndex, pivotIndex - 1)
Call QuickSort(arr, pivotIndex + 1, lastIndex)
End If
End Sub
Sub Swap(ByRef a As Variant, ByRef b As Variant)
Dim temp As Variant
temp = a
a = b
b = temp
End Sub
呼び出しは次のようにする
Dim arr As Variant
arr = Array(4, 2, 3, 1, 5)
Call QuickSort(arr, 0, UBound(arr))
ノートパソコンをUbuntuにしたので、Ubuntuでいろいろやってみる。ノートPCのほうが消費電力少なそうなので、常時運用にはこっちのほうがいいのか、なんて思っている。
簡単なスクリプトなら、Raspberry Piで十分だけども持ってないので。
# -*- coding: utf-8 -*- """ Spyder Editor This is a temporary script file. """ import urllib3 import datetime import os diff_day = 1 #if yester day, 1 root_url = "http://www.jma.go.jp/jp/radnowc/imgs/radar/" area_codes = ["000","201","202","203","204","205","206","207","208","209","210","211","212","213","214","215","216","217","218","219"] hour_list = ["00","01","02","03","04","05","06","07","08","09","10","11","12","13","14","15","16","17","18","19","20","21","22","23"] min_list = ["00","05","10","15","20","25","30","35","40","45","50","55"] def download(url, file_name): connection_pool = urllib3.PoolManager() resp = connection_pool.request('GET', url) f = open(file_name,'wb') f.write(resp.data) f.close resp.release_conn() now = datetime.date.today().strftime("%Y%m%d") yesterday = datetime.date.today() - datetime.timedelta(days=diff_day) print(yesterday.strftime("%Y%m%d")) i=0 download("http://www.jma.go.jp/jp/radnowc/imgs/radar/210/202008081425-00.png","testtest.png") for areas in area_codes: new_path = areas +"/" + yesterday.strftime("%m%d") if not os.path.isdir(new_path): print(new_path) os.makedirs(new_path) for hours in hour_list: for mins in min_list: i += 1; for areas in area_codes: url = root_url + areas +"/" + yesterday.strftime("%Y%m%d")+ hours + mins + "-00.png" print(url) print(areas +"/" + yesterday.strftime("%m%d") + "/{:04d}".format(i)) download(url, areas +"/" + yesterday.strftime("%m%d") + "/{:04d}".format(i)+".png")
上記のコードを作って保存しておく。とりあえず名前はradar.pyとしておく。
で、定期的に
python radar.py
を実行できればよい。ただ、上記で実行すると、ルートフォルダ以下に作ってしまうので、次のようにディレクトリを移動してから実施というようなプログラムが必要。次の内容をradar.bhとでもしておいて
#! /bin/bash
cd /home/username/radar
python3 radar.py
cd /home/username/radar
python3 radar.py
hoge
crontab -eで
1 0 * * * radar.bh
などと登録するとよい。Gnu nanoで編集するとき、^Xなどの表示があるが、何だろうとおもったら、^XはEsc Esc xの順番にキーボードを叩けばいいらしい。文字入力と制御を切り替えるエスケープがEsc Escなんだろう、きっと。
気象衛星の画像を定期的に集めるのにいい方法がないかと考え、pythonでコードを作り、一日一回動かす、みたいなことを考えた。
でかけたときでも勝手に回収できるようにlinuxでも動かせるようにしておくといいかも。
以下のコードで一日分のデータを集められる。日付が変わった頃合いにこれを動かせば前日分のデータを集めきれる。
# -*- coding: utf-8 -*- """ Created on Sat Aug 8 14:10:25 2020 @author: m """ root_url = "http://www.jma.go.jp/jp/radnowc/imgs/radar/" area_codes = ["000","201","202","203","204","205","206","207","208","209","210","211","212","213","214","215","216","217","218","219"] hour_list = ["00","01","02","03","04","05","06","07","08","09","10","11","12","13","14","15","16","17","18","19","20","21","22","23"] min_list = ["00","05","10","15","20","25","30","35","40","45","50","55"] import urllib3 import datetime import os def download(url, file_name): connection_pool = urllib3.PoolManager() resp = connection_pool.request('GET', url) f = open(file_name,'wb') f.write(resp.data) f.close resp.release_conn() now = datetime.date.today().strftime("%Y%m%d") yesterday = datetime.date.today() - datetime.timedelta(days=1) print(yesterday.strftime("%Y%m%d")) i=0 download("http://www.jma.go.jp/jp/radnowc/imgs/radar/210/202008081425-00.png","testtest.png") for areas in area_codes: new_path = areas +"\\" + yesterday.strftime("%m%d") if not os.path.isdir(new_path): print(new_path) os.mkdir(new_path) for hours in hour_list: for mins in min_list: i += 1; for areas in area_codes: url = root_url + areas +"/" + yesterday.strftime("%Y%m%d")+ hours + mins + "-00.png" print(url) print(areas +"\\" + yesterday.strftime("%m%d") + "\\{:04d}".format(i)) download(url, areas +"\\" + yesterday.strftime("%m%d") + "\\{:04d}".format(i)+".png")
Irvineという画像をダウンロードするプログラムがある。
気象レーダー画像をダウンロードするのに使っているが、ファイル名が規則的であるので、未来のファイル名も把握できる。
平時は忙しくてかき集めるのを忘れがち（連続で保存してたのに、１０枚途切れた、とか悲しい）なのだが、irvineというソフトに放り込んでおいて定期的にダウンロードしたらいいのではと思った。
まず、Irvineは下のようなソフトである。今回ダウンロードするのは気象庁の雨雲レーダー画像です。
ファイル名は、地域にもよりますが、日本全国の場合は
http://www.jma.go.jp/jp/radnowc/imgs/radar/000/202007220015-00.png
というように、ファイル名はyyyymmddhhmm-00.pngとなっている。時間は5分ごとである。このファイルリストを作成してIrvineに放り込む。
文字列は過去に作った文字列連結のJavaScriptで作っておく。
８月は31日あるので8928データが出力される。これをコピペでOK。
次からIrvineのスクリプト。
オプション設定からスクリプトを選ぶ。
もう作ってあるが、毎日ダウンロードというものがある。
ヘルプファイル（最近のWindowsはhlpファイルに対応してない！）を見ると、Start()というものが開始の制御のようである。
プログラムの中身は上記の通り。86,400,000msが一日なのだが、なんとなく半分の43,200,000msを周期にしている。
たぶん、これで一日に二回、一斉開始のコマンドとなるはずなので数日後に収穫を確認してみよう。
そして、気がついた。フォルダごとにタイマー設定が可能だった。スクリプトなしで。
たぶん、このタイマの開始タイマが機能するだろう。
測定データと近所の気象庁データとの比較をしてみた。
おおよそ600Pa程度ズレがあるが、傾向は同じである。ちなみに600Paずらしてみてもピッタリ重ならない。まあ、しょうがない。
値のズレが大きい（標高差でも2hPa程度のズレ）なので、調べると部屋を占めると24時間換気の影響で屋内陰圧になることがわかった。半ばわかっていたけど。
別の日の測定結果を見るとこのような形で気象庁の測定データに近い値担っている。せいぜい0.5hPa程度の差に収まっている。
マイコンでアナログデジタル変換（ADC）を行うとき、ローパスフィルタを入れたち時がある。でも部品が増えるの微妙ってときもある。
また、デジタルフィルタだと、特定の区間だけサンプリングしてフィルタするなどが容易である。
nをフィルタに使う整数、xを入力データとして、x1, x2, x3・・・・とあるものとします。y1=x1として、あとは上記の式に従って計算するとi番目の出力データyiが計算できます。
なお、整数演算する場合などは丸め誤差に注意が必要です。浮動小数点を問題なく扱えるマイコンなら、この限りではありません。
実際にこのフィルタを使った事例について。
関数は初期値30、途中から80になるというステップ関数に、±5の範囲でノイズを与えています(ただしエクセルのrand()関数なので分布が正規分布ではありません）。
関数だと、y=30+50*stepFunction(x) + (rand()-0.5)*10 といった感じ。
n=1は要するにフィルタなしです。上記の式にn=1を入れるとyi=xiとなります。
n=2～n=64まで比較するとアナログフィルタのように、応答が遅れ、終端部が平滑になることがわかります。
最後の100個分のデータを見ると次のよう担っています。
100データの標準偏差をとってみるとわかりますがnが大きいほどばらつきが少なくなります。
8bitマイコンではADCは10bitか8bitくらいなので、16bit変数を用意しておき、
z[i] = (n -1)y[i-1] + x[i] = ny[i-1] - y[i-1] + x[i] = z[i-1] - y[i-1] + z[i]
y[i] = z[i]/n
と、中間にz[i]というn倍大きな値の変数を入れておくと演算の負荷が多少すくなるなる（浮動小数点演算よりは計算がかんたん）。
たとえば2, 4, 8, 16, 32など2^nのサンプリングの場合はビットシフトで掛け算が済むので、すこし早く処理ができます。
ふと一日の気圧変化を記録に残そうと思ったときに、Arduinoには正確な時間計測の方法がなさそうだとなり、Processingで測定データを受信し、時刻を付与してファイルに記録、といようなことをしてみた。
Arduino側プログラム
MPL3115Aに設定して値を読み込む。今回は標高に換算したデータを出力するようにしています。
タイマー割り込みを使って0.25秒に一回、データをシリアル通信で吐き出しています。もっと遅くてもよかったと今では思ってる。
#include#include #define OverSample 128.0 #include "SparkFunMPL3115A2.h" //Create an instance of the object MPL3115A2 myPressure; LiquidCrystal lcd(12, 11, 5, 4, 3, 2); float Altitude = 0; float Ondo = 0; float Pressure = 0; float tempAltitude = 0; ISR(TIMER1_COMPA_vect) { digitalWrite(13, !digitalRead(13)); Serial.begin(9600); // Start serial for output Serial.print(Ondo,3); Serial.print(" "); Serial.println(Altitude,3); Serial.end(); } void setup() { pinMode(13, OUTPUT); TCCR1A = 0; TCCR1B = 0; TCCR1B |= (1 << WGM12) | (1 << CS12); //CTCmode //prescaler to 256 OCR1A = 156250-1; TIMSK1 |= (1 << OCIE1A); lcd.clear(); lcd.begin(16,2); lcd.setCursor(0,0); lcd.print("Hello! Wait!"); lcd.setCursor(0,1); lcd.print("OverSample:"); lcd.print(OverSample); Wire.begin(); // Join i2c bus // delay(100); // Configure the sensor myPressure.setModeAltimeter(); // Measure altitude above sea level in meters //myPressure.setModeBarometer(); // Measure pressure in Pascals from 20 to 110 kPa myPressure.setOversampleRate(7); // Set Oversample to the recommended 128 myPressure.enableEventFlags(); // Enable all three pressure and temp event flags } void loop(){ myPressure.begin(); // Get sensor online for(int i = 0; i< 10; i++){ Altitude = myPressure.readAltitude(); // Pressure = myPressure.readPressure(); Ondo = myPressure.readTemp(); } while(1){ Altitude = Altitude * (OverSample-1)/OverSample + myPressure.readAltitude()/OverSample; // Pressure = Pressure * (OverSample-1)/OverSample + myPressure.readPressure()/OverSample; Ondo = Ondo* (OverSample-1)/OverSample + myPressure.readTemp()/OverSample; lcd.clear(); lcd.setCursor(0,0); lcd.print(Ondo, 2); lcd.print("C"); lcd.setCursor(7,0); lcd.print(Altitude, 2); lcd.print("m"); } }
Processing側プログラム
import processing.serial.*; PrintWriter file; Serial myPort; // set Serial port to String inString; // string to read from serial port int lf = 10; // ASCII \n void setup() { file = createWriter("balomatic_tester"+nf(year(),4)+nf(month(),2)+nf(day(),2)+"_"+nf(hour(),2)+nf(minute(),2)+".log"); size(200,200); printArray(Serial.list()); myPort = new Serial(this,Serial.list()[1], 9600); myPort.bufferUntil(lf); } void draw() { } void serialEvent(Serial p) { inString = p.readString(); String now; now = nf(year(),4)+"/"+nf(month(),2)+"/"+nf(day(),2)+" "+nf(hour(),2)+":"+nf(minute(),2)+":"+nf(second(),2)+" "; print(now); print(inString); file.print(now + inString); } void keyPressed(){ println(key); if(key == 'f'){ file.flush(); } else if(key == 'q'){ file.flush(); file.close(); exit(); } }
void keyPressed()でやっているのは、fキーを押したときにはファイルをフラッシュ（ディスクに書き出す）を実施、qをクリックしたら実行中止、というようなことになっています。
これでProcessing側でデータを受信して、ファイルに吐き出していく。コンソールには次のようにログが出てきます。
アクセス
トータル
