C++/STLでDFT、IDFTしてみる。

C++/STLでモダンに複素DFT、IDFTを実装。なんか、もやもやしていたけどすっきりした。やってみたら簡単だった。


#define _CRTDBG_MAP_ALLOC
#include <stdlib.h>
#include <crtdbg.h>

#include <iostream>
#include <complex>
#include <cmath>	// sqrt, sin, cos
#include <iomanip>	// setprecision

const double pi = 3.1415926535897931f;

const int length = 32;

// initialize array
template <class T> 
void inita(std::complex<T> *data, int length) {
	for(int i=0; i<length; i++) {
		data[i] = std::complex<T>(0.5f+std::sin(2*pi*i/length+4), std::cos(4*2*pi*i/length));
		//data[i] = std::complex<double>(2, 3);
	}
}

// display array
template <class T>
void dispa(std::complex<T> *data, int length) {
	for(int i=0; i<length; i++) {
		std::wcout << i << L":" << std::fixed << std::setprecision(5) << data[i] << L" " << abs(data[i]) << std::endl;
	}
}

// deep copy
template <class T>
void copya(std::complex<T> *data, int length, std::complex<T> *dest) {
	for(int i=0; i<length; i++) {
		dest[i] = data[i];
	}
}

// DFT
template <class T>
void worka(std::complex<T> *data, int length, std::complex<T> *res) {
	for(int j=0; j<length; j++) {
		res[j] = std::complex<T>(0.0f, 0.0f);
		for(int i=0; i<length; i++) {
			res[j] += exp(std::complex<T>(0,  2*pi*i*j/length)) * data[i];
		}
	}
}

// IDFT
template <class T>
void workb(std::complex<T> *data, int length, std::complex<T> *res) {
	for(int j=0; j<length; j++) {
		res[j] = std::complex<T>(0.0f, 0.0f);
		for(int i=0; i<length; i++) {
			res[j] += exp(std::complex<T>(0, -2*pi*i*j/length)) * data[i];
		}
		res[j] /= length;
	}
}

void main(void) {
	_CrtSetReportMode(_CRT_WARN,  _CRTDBG_MODE_FILE);
	_CrtSetReportMode(_CRT_ERROR, _CRTDBG_MODE_FILE | _CRTDBG_MODE_DEBUG);
	_CrtSetReportFile(_CRT_WARN,  _CRTDBG_FILE_STDERR);
	_CrtSetReportFile(_CRT_ERROR, _CRTDBG_FILE_STDERR);
	//#define d_new new(_NORMAL_BLOCK, __FILE__, __LINE__)

	using atype = double;
	std::complex<atype> *data, *res;
	data = new std::complex<atype>[length];
	res  = new std::complex<atype>[length];
	
	std::wcout << L"-- orig --" << std::endl;
	inita(data, length);
	dispa(data, length);
	
	std::wcout << L"-- DFT --" << std::endl;
	worka(data, length, res);
	copya(res, length, data);
	dispa(data, length);
	
	std::wcout << L"-- IDFT --" << std::endl;
	workb(data, length, res);
	dispa(res, length);

	delete [] data;
	delete [] res;
	
	_CrtDumpMemoryLeaks();

	return;
}

ビルドはVC++で cl /EHa complex.cpp

十数年前に学校で履修してわけがわからないまま単位取得したテイラー展開フーリエ展開ラプラス変換、それと線形代数の行列・逆行列はすべて関連しているんだと、なぜか未明に気づいてはっとした。数学が得意な人からしたら当たり前じゃんてなるだろうが、数学が嫌いな人からしたらそれぞれは離散的で結びつかない。こういうことはもっと早く教えてほしかった。ウィキペディアの説明ではとても理解できない。



忘れないうちにコードのメモ。

main()内の_CrtSetReportModeあたりはCランタイム組み込みのメモリリーク検出デバッグ機能の有効化。最後の _CrtDumpMemoryLeaks();を実行するときに開放されていないメモリがあれば警告してくる。#define d_newを有効にしてヒープのnewをd_newに変更すればどこでメモリリークしたか表示する。

using atype = のところで全体の精度を指定する。どうせfloatかdoubleしか指定しないが、一応テンプレート対応てことで。Old C Styleのtypedef ~ と同じ。length = の指定で長さは任意に変更できる。FFTではないのでn^2長さでなくても構わない。

inita()で初期データを作成する。サンプルはY=0.5+SIN(2πt+4) + COS(4*2πt) i のオフセットのある実数部正弦波+虚数部余弦波のツートーン。

worka()、workb()がDFT、IDFTのメイン処理。同じデータ長の実数部と虚数部で90度位相差がある二つの正弦波を0倍(DC)~N-1倍の振動数まで順に掛け合わせる。DFTとIDFTは掛け合わせる振動の向きが正負逆になる。オイラーの公式で実数と虚数に分離してSIN、COSの組み合わせで計算すると逆行列になる。

プログラミング書籍なら翔泳社の通販『SEshop』

いちばん簡単でいちばん詳しいDX-SR9Mのデジタルモード付きな令和最新の局免申請のし方。

アルインコのHFオールモード名機DX-SR9M/J

高めのハンディ機くらいの値段なのに、ドットマトリクス液晶やDSPを内蔵せず、最新のアナログ技術を結集して作られた、最後にして最高の大手メーカー製アナログ機といえよう。起動時のアルインコ音が無いのもアルインコっぽさが無くて良い。各モード用の標準・ナロー ムラタフィルター、R8ゆずりのLF~のゼネカバ受信、コリンズフィルターを任意に搭載可能、TCXOで1ppm、そしてSQL・RIT・IF-shiftを自由に論理位置変更できる謎機能搭載。ごちゃんねらーですら文句のつけようがほぼ無い仕上がり。まさに名機。

名機なのに技適無しで申請が面倒。名機なのにネットを探しても開局申請したという報告がほぼ無い。6エリアのOM局長様がTSSでデジタルモードを保障認定したことを報告されているが取説通り・メーカーの系統図で、ということになっている。アマチュアのデジタルモードは去年簡素化されたので取説とメーカーの開示している情報通りだと簡略化前の内容になっている。

12年以上前に当局のデジタルモード申請をしたときはTSSやら総通から問い合わせが来て回答書を作ったりして非常に面倒だった。さすがにデジタルモードが浸透した今はそんな面倒な事にはならないだろうけど、今回もなるべく楽したいし、附属装置を少しいじったくらいで再度申請とか必要無いように申請したい。

ということでDX-SR9Mを使ってデジタルモード有りで、できる限り簡単な内容で保証認定と免許申請をしてみる。今回は初めてJARDの保証認定を使ったがTSSでも同様だと思う。ちなみにアルインコが商品ページでTSSではなくJARDを推しているのはAlincoがJAIA加盟メーカーだからと推測した。JAIAの関連団体リンクにあるアマチュア無線関連団体の順列ではJARDが先頭にあってTSSが一番最後にあるからJAIAの2社に対する扱いが読み取れる。アルインコの商品ページでJARD推しになってるのは、TSSよりもJARDを使ってくれというアルインコの隠れたメッセージだろう。(並びは50音順だろうと言われればそうかもしれないが、JARDとJAIAは関係性が深い。)

局免申請の「事項書及び工事設計書」
13 電波の型式並びに希望する周波数及び空中線電力

3アマさんは10M、14M無しで。

15 備考 簡素化デジタルモード手続きのキーワード「デジタルモードのため附属装置(PC)を接続」を書いておく。(実際の申請では “PC等” って書いた。ICレコーダーとか繋ぎたいし。)
16 工事設計書 工事設計書の中は後述。
添付書類 JARD/TSSから送られてきた保証認定の書類を追加。

16 工事設計書の工事設計書編集から工事設計情報入力を行う。基本的には取説から拾って次の項目を埋めるが、取説からは変調方式の入れ方が分からない。画像参照。
・発射可能な電波の型式及び周波数の範囲:事項書の一括表記を超えない範囲で。
・変調方式
・終段管
・定格出力
・添付書類:後述

電波の形式を入れるのが面倒。
3アマさんは10M、14M無しで。

添付書類には送信機系統図を追加する。画像ではPDFを添付しているがJPEG二枚でも良いはず。
1ページ目は取説をスキャンしてDX-SR9M(50W機)の内容だけにしたもの。取説はDX-SR9J(100W機)と共用になっている。
2ページ目は附属装置の接続を書いたもの。これが今回のポイント。メーカーの説明のようにモードごとの細かい諸元は書かない。

スキャンした画像を少し編集する。
コリンズフィルター装着の場合はフィルターのところも編集。

取説やメーカーの説明ページにある送信器系統図は簡素化されたデジタルモードの手続きから照らし合わせると冗長すぎるし、将来新モードを追加する場合は申請内容との不整合で疑問が生じる。

送信器系統図の2がやたら簡単になったが保証認定、免許申請は通る。



DX-SR9Mは最初から50Wに設定された送信機だから3アマでも問題なく保証認定されるが、DX-SR9Jを50Wに設定した送信機を3アマが申請する場合はメーカーかJAIA加盟店(販売店)が作成した「空中線電力50W固定措置に関する証明書」を添付する必要があるかもしれない。1アマ・2アマの場合は改造箇所の写真とパワー計で各バンドの出力を測定している写真を添付すればOK。移動する局の申請で面倒に巻き込まれたくなければ1アマでも最初からDX-SR9Mを買う。

名機なのに高めのハンディ機くらいの値段しかしないから誰もが買いやすいのに10W機の設定が無く、残念ながら4アマはこの名機を使うことはできない。仮に送信ができる状態のDX-SR9M/Jをアンテナを繋ぐと不法無線局の開設に対する罰則として「1年以下の懲役または100万円以下の罰金」が科せらるほどの犯罪行為。3アマか2(メール)アマにステップアップしてこの名機を堪能しよう。



10月_3日 保証願提出
10月21日 保証完了
10月21日 変更申請提出
12月_1日 変更申請審査終了
12月_2日 返送用封筒送付
12月_8日 免許状到着

今回は一回の質疑も手戻りも無く完了。
かかった日数67日。長かった。

スタンド充電器の中を確認。

ICOMのハンディトランシーバー用の急速充電器、スタンド型だから置いて使うときなど便利なのに基本別売でそこそこ高価だからオプション商法などと揶揄される。

現在唯一のハンディ機ICOM IC-T70もそのパターンでこの機種専用のBC-191が入手できなくなる前に買った。

まあ、ごく普通。

アイコム製の機器はこのあたりのカッチリ感がアルインコ製より良い感じ。アルインコのDJ-X7なんかはスタンドチャージャーが附属してるけどなんかグラグラしてて接触不良気味になったりする。

BC-191は本体ときちっと勘合して気持ちがいい。でも値段が高いのに中身が接点だけとか有り得ないし、とりあえず分解してみた。

回路のトレースはしたくないくらいの集積度。ルネサスの8ビットマイコンμPD789112、富士通のレギュレータIC MB3759、JRCのオペアンプNJM2904 なんかの国内メーカーの集積回路が載っている。基板も多層基板で受動パーツ類もしっかりした印象。アマチュア無線機用としてはもったいないくらいの作り。

アイコムのメイド・イン・ジャパンへのこだわりを感じた。

摺動式変圧器の中を観察する。

東芝の商品名スライダックで一般に通じる摺動式変圧器、よく見るものは円筒形のケースに入って調整ツマミが上面にあるタイプ。

手元にあるのはそれとは違って角型のケースに入って電流計、電圧計、スイッチ、ヒューズが内蔵されたもの。メンテナンスのために中を開けてみる。

山菱電機では「ボルトスライダー」と呼称する。東芝はスライダックの後継として山菱のボルトスライダーを推奨しているが、古くからある技術で作られるものなので他のメーカーでも品質に大差は無い。スペックは入力100V、出力0~130V 15Aでごく普通。トランス本体は横に固定されている。重量があるのでゴムのスペーサーで下、左右を支えている。電流計、電圧計はともに信頼の横河電機製でCLASS 2.5のJIS規格品。

スライドトランスの製品で電流計とヒューズを内蔵したものは結構少なくて、電圧はダイヤルに目盛りが付いていればだいたい分かるものの、電流は電流計やクランプメータを使わないと全く分からない。それに電流ヒューズを内蔵していないものも多くて安全に使おうと思うと木板に円筒形のスライドトランスと電流計、電圧計とサーキットブレーカかヒューズを固定した治具を作ることになる。これが結構かさばって取り回しが良くないし上に物を置けないから場所を取るしでスマートじゃない。その点、既製品で角型のスライドトランスは高価だけれど便利。

電圧降下でトランスが必要になる状況なんてよっぽど大電流を長い引き回しで使うのでなければ、商用電源が安定している現在の一般家庭では必要無いが、ちょっとした電気の実験をしたいときに昔ながらのスライドトランスは一台あると助かる。

参考: 直流安定化電源を分解してみた。
https://mzex.wordpress.com/2019/09/07/14838/

『顔と体のための個人電気分解』を分解してみた2。

家庭用電気針脱毛器を分解したら、予想していたよりおもしろい回路になっていた。

前回記事:
『顔と体のための個人電気分解』を分解してみた。

Clean+Easy Home Electrolysis Kit Machine Hair Removerの回路をトレース。一部は意図的に削除。

シンプルに見えるけど動作を追うと結構よくできている。汎用ロジックIC 1個でこんなにいろいろ実現しているのはすごい。毛抜きの内部でこんなに感動できるなんて思ってもみなかった。

・基本的な動作として、待機状態では針から皮膚に流れる電流は微量で痛み無しに針を毛穴に挿入できる。挿入した針が毛根に到達して電気抵抗が下がり針に流れる電流が60μAを超えると動作状態に入り電気分解に必要な0.4~0.8mA程度の電流が針から約8秒間供給される。その後待機状態に戻る。

・全体の電源は006P電池に逆接続防止のD1を経由したうえで供給している。(BATT 9[V] – VD1 0.6[V] で Vdd=8.4[V])
・微弱電流により抵抗の変化を監視する待機状態と、電気分解を行う電流を流す動作状態がある。
・待機状態から動作状態の遷移はR1とスタイラスが当たっている抵抗値の比率でヒステリシス付きNANDゲート1/4の入力がH→Lのスレッショルドを下回った時に発生する。
・動作状態から待機状態への遷移は、C2とR4で構成される時定数回路で決定される。約8秒。
・スタイラスのプレートP1が+で針はGNDにつながっている。P1は100kΩを経由しており、待機状態では動作状態に遷移する直前で60μA程度しか流れない。(NANDのH→Lのスレッショルド電圧VH→LがVddの33%だと仮定して Vdd / (100k+50k)[Ω] = 56[μA]) この程度だと皮膚は刺激を感じない。
・本体のプレートはスタイラスのプレートと完全に同電位ではなく、R8を経由しており動作状態でのショートによる針先の破損を保護する。スタイラスのプレートと針は構造的にショートしないのでR8の制限が必要ない。
・スタイラスの先端と本体のプレートを直接接触させることで動作状態に遷移してブザーが鳴るためこれを動作確認として利用できる。
・スタイラスを毛根に挿入してR1の電圧降下が発生するとR2の電圧も下がり、NAND 1/4の閾値を下回ると出力がHになり時定数回路を構成するC2の充電を開始、C2の充電が行われる時間がNAND 2/4のワンショット動作の時間となる。
・NAND 2/4のワンショット動作が始まると、NAND 4/4の出力がHになり、この電圧が操作者が設定したVRを経由してスタイラスのプレートP1に電気分解可能な電流が供給される(動作状態)。このとき最大電流はImax[A] = (Vdd – VD3)[V] / (VR+人体抵抗)[Ω] となり、VRはダイヤルの位置により10.4kΩ~1.5kΩ程度の範囲で可変される。仮に人体の抵抗を8kΩ程度としたら、VRの位置により電流はおよそ0.4mA~0.8mAの範囲となる。
・最大電流はNANDゲートの出力電流の制約によりVRのショートモードでも10mA程度に制限される。(そこまでは想定していない?)
・ダイヤルは半固定抵抗に直結されているだけなので0の位置でも動作状態に電流は流れる。
・C4の充電が完了するとNAND 2/4の出力が再びHになり、NAND 4/4がオフになってスタイラスへの動作状態の電流が遮断される。C1とR6による遅延回路でNAND 2/4と4/4の連動はわずかに遅れるためNAND 2/4のワンショット動作終了時はNAND 1/4を確実に待機状態へ遷移できる。
・待機状態に遷移後も電気抵抗が低い状態であれば再び動作状態を開始する。
・C2充電中はNAND 3/4で構成される発振回路で圧電ブザーを駆動する。
・待機状態でもわずかに電源を消費するため未使用状態でもダイヤルの位置に関係なく電池は無くなっていく。

昨日紹介したリンク先の情報から、抵抗検出つきの電気分解法回路は特許取得されており、回路は異なるが公開されている。この特許について現在は期限切れ。今回購入したPersonal Electrolysisにパテントの記載は無い。現物の回路は特許文書とはだいぶ異なっていて、この特許文書の回路は挿入時適正な電流値になったときから一定時間ブザーを鳴らす機能だけとなっている。また電流の検出方法も違う。

Electrolysis hair removal apparatus – Google Patents
https://patents.google.com/patent/US4216775A/en


毛抜き無しバージョンの回路も推定してみた。

分解の写真や前回紹介したリンクから推定する限り、これだけだと思う。特許回路は実装されていない。回路としては自作バージョンと同じ。ボリュームは10kΩらしい。毛抜き有りバージョンと比較して手抜き感が半端ないコストダウンの跡がうかがえる。工夫があるとしたら、ケースの物理的な構造によりボリュームの動く範囲が制限されて、操作を間違えても極端な電流がかからないようになっているくらいなのと、使っていないときは全く電源を消耗させないことくらいか。毛根に挿入するまで正極側を触らないようにしないと、毛根に達する前に電気分解に必要な電流が供給されてしまって痛いし操作が面倒。今は毛抜き無しバージョンが売られてないから間違わないけど、もし選ぶなら毛抜きの差だけではないから毛抜き有りバージョンにしよう。

※毛抜き有りバージョン、毛抜き無しバージョンについて。
販売元のAmerican International Industries社はブランド名としてclean+easyで毛抜き有り版だけが出ているが、以前はOne Touchブランドでエレクトロリシス製品を出していて、毛抜き有りがDelux版、毛抜き無しが無印版で出ていた。One TouchブランドのDelux版はclean+easyの現行品と全く同じ構成。毛抜きの有り無しの違いだけかと思ったら回路が全く違ってることが判明した。

ちなみに付属の毛抜きは極めて使いにくい。外国のレビュー動画でも「毛抜きはたくさんあるからもう要らないヮ」みたいなことを言っていた。