超音波風速計

　ずーっと前からやりたくて、ちょこちょことやってたんだけど、また久しぶりに超音波風速計ネタ。

　試しにSTBee Miniで40kHzのPWMを作ってADCでサンプリング。PWMは位相が180度ずれた2chを出力し、6.6Vppを超音波スピーカーに入力。受信は超音波マイクからADCにダイレクトで接続。

　1枚目の画像は4回のサンプリングの波形、2枚目は、1枚目の最初の方を拡大した図。PWM生成のバグというか不具合というかで位相が逆になってしまっているが、正確に180度ずれているだけで、余分なずれは見られない。
　3枚目は、先の4計測と、マイク側とスピーカー側から息を吹いた、計6回分の波形。若干のオフセットが有るが、無風状態では完全に位相差がゼロで、息を吹いたときは位相差として計測できているのがわかる。
　縦軸は電圧(出力値*3.3・4096で計算)、横軸は時間で単位はマイクロ秒。


　先に、「マイクからADCにダイレクトで」と書いたが、本来はオペアンプで増幅してからADCに入れようと思っていた。しかし、手持ちのオペアンプでは帯域が足りないらしく、アッテネータにしかならなかった。しょうがないから試しにマイコンに直結したら、0.5Vppくらいでサンプリングできてしまった。
　送信もパスコンでACカップリングしただけの直結だし、受信もマイコンに直結だから、今回は能動的なアナログ回路は一切無い。実環境で使えるかはわからないが、風速数m/s程度ならデジタル回路オンリーでも意外となんとかなりそうな気がする。
　正直あんまり複雑な基板を作る気力がないので、てきとーに配線するだけでも動きそうな、アナログ回路レスな基板というのは結構魅力的。希望が出てきた。


　現在使っているハードウェアリソース：TIM2(PWMタイミング)、TIM3(ADCタイミング)、DMA1-1(ADC転送)、DMA1-2(PWM転送)、ADC1。
　TIMとDMAとGPIOでPWMを半ハードウェア的に作っているので、PWMチャンネル数は16までいける。スピーカー1個に2chを使うから最大8個まで。アナログ入力はGPIOAだけだと8chまで可能。GPIOは大量に使うが、SP8ch、MIC8chまでなら現在のリソースだけでいける。ただ超音波回りだけでGPIOを24本も使うと、他に自由に使えるのはUSB回りを除いた5本くらいしか残らない。まぁUSART1は競合しないから、もう一個マイコンを載せて表示とかはそっちに任せるのもありだが。スピーカーからIOを2本もらってくればCANバスが使えるので、長距離を引き回して表示するならそっちのほうが良いかな。


　とりあえず今のままでも簡易的な1次元の風速計は作れると思うので、それが動いたらさらに軸数を増やしてみる。単純な信号処理のテストならPCで動くし、いちいちマイコンに転送したりしなくていいので楽。

FREEDOM

　CelestrakのTLEからFREEDOMが消えました。

　とりあえず例によって地心からの高度グラフ。

　黒矢印がついてるのが今日(2017/02/08)です。TLEから消える前の数日分の軌道要素がCelestrakに無いため、予測値を書いています。
　FREEDOMの放出はこのグラフで42751のあたりです。おおよそ25日程度で再突入したとみられます。他のキューブサットと比較して、かなり短い期間で再突入したといえます。

STM32F1でADCをTIMで駆動してDMAを駆動する

　やりたいこと：ADCで連続変換、タイミングはTIMで作ってDMA転送
　つかうもの：STM32F1

　今回はTIM4でタイミングを作り、ADC1を駆動、PB0のchを変換してDMA転送を行う。

　ADC1をTIM4で駆動する場合、TIM4 OC4のエッジで変換が開始されるが、注意点として、OCのOutputStateはEnableにしなければいけない。
　例えばTIMをトリガにDMA転送を行う場合、内部的にOCのエッジがあればいいので、OutputStateはDisableでも問題ない。しかし、ADCのトリガとして使うにはOutputStateはEnableにしなければ動かなかった。
　たとえOutputStateをEnableにしても、GPIOをAFで初期化しなければOCが出力されることはない。GPIOをソフトウェアで入出力に使う分にも問題はないが、ペリフェラルで使うことはできなくなる。TIM4 OC4はPB9に接続されており、このピンはI2C1 SDAのリマップと、CAN TXのリマップに接続されている。特にCANは標準がUSBと同居しており、STBeeMiniで使う際はPB8/PB9にリマップする必要があるから、TIM4 OC4とは競合してしまう。もしかしたらTIM4 OC4をリマップしてPD15に移動すれば問題ないかもしれないが、そこまでは確認していない。


　TIMでADCを駆動してDMA転送すれば、バッファが埋まるまではCPUからはノータッチで処理できるので、他の処理に専念することができる。YES DMA! NO タッチ!

高速じゃないけどちょっと高速なフーリエ変換

　正数演算で離散フーリエ変換してみた。

void DFT(const int* f, int N, int *result) {
    int n, nEnd, k;
    nEnd = N /2;

    for (n = 0; n < nEnd; n++) {
        int re = 0;
        int q = n * 1000;
        for (k = 0; k < N; k++) {
            int p = q * k / N;
            re += f[k] * CosTable[p % 1000];
        }
        result[n] = re;
    }
}

　他に1キロバイト分の余弦波テーブル(int8_t CosTable[1000])が必用。テーブルは0.01フェーズステップで、エクセルとかで作っておけばok。ROMを節約したい場合は0.25フェーズ分で250バイトでも良いけど、実行時間とのトレードオフになる。
　虚数は計算しないので、離散フーリエ変換の半分のみだが、それでもO^2/2の計算量なので結構時間がかかる。例えば高速なフーリエ変換と比較して、128ポイントならおよそ9.2倍、512ポイントなら30倍、1024ポイントなら52倍くらいの差がある。とはいえ、72MHzのARMで128ポイント計算するのに必用なのは5msec未満程度だったので、ポイント数が少ない場合はFFTのややこしいアルゴリズム通すより楽だと思う。あとポイント数が自由に選べるというのも地味に嬉しいかも？ま、FFTのほうが明らかに早いんですけどね。

比較的高いところまでの気圧とか

　なんとなく成層圏とかあのあたりの気圧が欲しい気分だったので計算してみた。

　縦軸が高度(km)、横軸は上が気圧(hPa)、下が気温(degC)で、赤が気圧、青が気温。

　ググって出て来るグラフに近いので、そんなに大きくは外してないと思う。本来はグラフは85kmあたりで終わりなのだが、このグラフではそれ以上も続いてる。

　値は国際標準大気の米wikipediaページを使った
　International Standard Atmosphere - Wikipedia
　ちなみに日本語ページと英語ページでは微妙に値が違う。今回はだいたいそれっぽい値だったら良いというスタンスなので気にしない。

　参考までに、C#のソースコード。気圧の計算は以下のページの式をそのまま使った。
　目的地の気温と気圧を計算 - 高精度計算サイト
　レイヤの計算がかなり冗長なのであまりにも大量に呼ぶ場合は最適化したほうが良いかも。

static void CalcSeaLevelPresAndTemp(double TargetAltMeter, out double TargetTempDeg, out double TargetPresHPa)
{
    int[] BaseAlt = { -610, 11000, 20000, 32000, 47000, 51000, 71000, 84852 };
    float[] BaseTemp = { +19.0f, -56.5f, -56.5f, -44.5f, -2.5f, -2.5f, -58.5f, -86.28f };
    float[] BasePres = { 1089.0f, 226.32f, 54.749f, 8.6802f, 1.1091f, 0.66939f, 0.39564f, 0.003734f };
    float[] TempLapseRate = { -6.5f, +0.0f, +1.0f, +2.8f, +0.0f, -2.8f, -2.0f, +0.0f };

    int Layer =
        TargetAltMeter >= BaseAlt[7] ? 7 :
        TargetAltMeter >= BaseAlt[6] ? 6 :
        TargetAltMeter >= BaseAlt[5] ? 5 :
        TargetAltMeter >= BaseAlt[4] ? 4 :
        TargetAltMeter >= BaseAlt[3] ? 3 :
        TargetAltMeter >= BaseAlt[2] ? 2 :
        TargetAltMeter >= BaseAlt[1] ? 1 : 0;

    double deltaAlt = TargetAltMeter - BaseAlt[Layer];
    double T = deltaAlt / 1000 * TempLapseRate[Layer] + BaseTemp[Layer];
    double P = BasePres[Layer] * Math.Pow(1 - ((0.0065 * deltaAlt) / (T + 0.0065 * deltaAlt + 273.15)), 5.257);

    TargetTempDeg = T;
    TargetPresHPa = P;
}

　うぃきぺのグラフだと高度85kmを上回ったあたりで切れてるので、レイヤ>=7で例外を投げるとかした方がいいかも。


　ちなみに理科年表のグラフによると、密度は高度100kmあたりまで対数でほぼ線形の減少となり、気温は100kmを超えたあたりで一気に加熱されている。高度400km以上は1000degKで変化せず、という感じらしい。宇宙の「気温」はだいたい摂氏700度くらいということなので、めっちゃ熱いはずなのだが、宇宙はかなり寒い印象があると思う。そもそも空気密度が限りなく低いので、熱伝導で加熱されることはなく、放射でのみ加熱・冷却されるから、気温の影響は無視できる程度なんであろう。


　ということで、なんとなくそれっぽい気圧とか高度を求める話でした。

スペクトラム拡散

　なんとなくスペクトラム拡散をやってみたくなったのでやってみた。

　上から準に無変換連続波、周波数ホッピング、直接拡散、のスペクトル画面。横軸が周波数、縦軸が時間。

　青が無変換連続波、赤が直接拡散のスペクトル。横軸が周波数なのは同じだが、縦軸は対数になってる点に注意。デシベルじゃなくてベル。

　無変換と直接拡散では信号強度(絶対値)に4倍ほどの差がある。ノイズフロアに埋もれるというほどじゃないけど、いろいろノイズの有る空間ならもうちょっと見分けづらくなると思う。

　無変調はスペクトルを見るまでもなくAM復調でスイープしていけば、信号が出ていれば一発でわかってしまう。
　周波数ホッピングでは連続で受信することはできないが、スペクトルを見てればホッピングの頻度や拡散がどれくらいかというのは見てわかる。最近だとソフトウェア無線で広帯域を一気に受信できるから、送信者が1なら追尾できそうな気もする。
　直接拡散では薄っすらと信号があることはわかるが、環境や符号によっては送信されてる事自体がわからなくなってしまう。

　昔のレーダーは無変換だったので、何らかの手段で周波数さえわかれば、その周波数にAMラジオをあわせておけばレーダーが照射されたことを把握できるし、逆位相を返すようなトランシーバーがあれば無効化できるし、タイミングを合わせて送信すれば居もしない影を見せつけることができるし、そんなことしなくたって強い電波を当ててやればレーダースクリーンを塗りつぶすことができる。レーダーが使われた初期の頃は電子部品も性能や信頼性の問題があったが、最近ではどんどん高性能化され、色々と都合が悪くなってきた。
　周波数拡散を行うことにより、拡散パターンが知られない限りはレーダーや通信を無効化される危険性が少ない。とはいえ、昔は「周波数を知られない限りはレーダーや通信を無効化される危険性が少ない」と言われていたはずだから、情報量を増やして相手が把握するのを遅らせる、程度にしかならない。数年・数十年も経てばそういうのに対抗する手段もいろいろ出てくるだろうし、さらにそれに対抗する手段も出てくるはず。
　直接拡散の場合は特定のスペクトルに集中するエネルギーが少ないから、適切な拡散符号で復調しない限りは送信されていることすらもわからない。戦闘機等のレーダー警報は強いエネルギーに対して警告を発するだけだから、周波数拡散を行ったレーダー波には反応できない。帯レーダーミサイルも同様で、強いエネルギーに突っ込んでいくだけだから、直接拡散とか使われるとつらそう。
　

　とりあえず今回は拡散だけだったので、次は復調とかもやりたい。とはいえプログラムの中で波形いじるだけなのでなぁ。
　個人で電波で遊ぶのはかなり大変だし、かといってスペクトラム拡散とかで遊ぶとなるとどういう遊び方になるんだろうか。可聴帯域かなぁ。うるさそう。ソリトンレーダーでも作るか。ソリトンレーダーって音波だと思ってたんだけど、調べてみると電波っぽい。そりゃそうだよなぁ。MGS4だと音って言ってた気がするんだが。

ふりーだむ

　最近ネタがないので例によってTLEウォッチです。

　ISSとSTARS-C、それからFREEDOMの地心からの距離をプロットしたグラフです。衛星番号41930がFREEDOMのはずですが、Celestrakではまた入れ替わってます。

　地心からの距離なので、グラフ下端の6371kmあたりが地表となり、上端の6871kmあたりが高度500kmとなります。100km/div(25km/div)です。
　横軸は1900年1月1日からの経過時間を示しています。1day/divです。

　値が上下に振動しているのは、ISSと言えども完全な真円の軌道ではないためです。今のところISSは焦点が中心から5kmくらいズレた楕円軌道で、±5km(10km pp)の振幅を示しています。
　本来、人工衛星を緯度経度高度で表した場合、この振幅は2つの正弦波を重ねた波形となります。これは地球が極方向につぶれた楕円体であるため、南北方向の移動時に直下の地面の高さが変わるため、相対的に人工衛星の高度の変化として出てくるわけですが、見づらいので今回は地心からの距離をプロットしています。

　ISSは400kmを少し上回ったあたりにいます。FREEDOM放出直後もだいたい同じくらいの高度ですが、最近では遠地点が350kmを下回り、近地点は325kmを下回っています。傾きとしては1日に10kmくらいのペースに見えるので、あと数日の間に近地点が300kmを下回り、直後には遠地点も300kmを下回りそうな感じです。

　以前、TRICOM-1の軌道を予測した時に、「まず遠地点が下がって円軌道に近くなり、やがて全体的に高度が下がる」と予想しましたが、このグラフを見る限りでは、その予想は外れています。
　FREEDOM放出直後の離心率はおよそ0.001ほどでしたが、最近では0.0017ほどとなり、少しずつ離心率が大きくなっています。上記の予想をした理由としては、遠地点よりも近地点のほうが大気が濃く、近地点のほうがより強く減速され、近地点での減速は遠地点の降下を引き起こす、という予想からでした。しかし高度400kmで膜を広げたような状態ではあまり関係ないようです(あるいはそもそもこの予想自体が間違っているか)。



　そろそろ違うネタも探してこないとなぁ。