<< LTSpiceによる歪み率の算出 | top | デイヴィッド・J・ハンド「偶然の統計学読了」 >>

スポンサーサイト

一定期間更新がないため広告を表示しています

スポンサードリンク | - | | - | - |

歪み率算出プログラム

先週から自転車通勤を始めた。
車が一台しかないので、昼間妻が車で出られるように。
あと健康のため。
車の運転が好きでないこともある。

会社まで3kmほどなのでそれほど大変ではない。
今の時期は気温もちょうど良く、快適な通勤だ。
渋滞も関係ないし。
冬になって自転車通勤する元気があるかが心配だ。

しばらく自転車に乗ってなかったから忘れてたけど、近所のコンビニとか行くのにめちゃめちゃ便利だね。
わざわざ新しい自転車を購入したので頑張って通勤します。

Fig. My自転車。かごがない方がかっこいいけど実用重視。下は砂利しきかけの防草シート。

///////////////////////////////
前回予告した通り、歪み率を算出するプログラムを紹介する。
時刻と電圧が記述されたファイルを入力として歪み率を計算する。
フーリエ変換はFFTWライブラリを使用した。シミュレーション開始直後は値が不安定なことがあるので最初の方は読み捨てている。
 
#include <stdio.h>
#include <stdlib.h>
#include <fftw3.h>
#include <math.h>
#define PI 3.14159265358979323846

int main(int argc , char **argv)
{
	int Nol = 0;//入力データ数
	double f = atof(argv[2]);//入力信号の周波数
	double S = atof(argv[3]);//サンプリング周波数 データから算出してもいいがとりあえず手入力
	FILE *fp = fopen(argv[1] , "r");
	double t , v;
	while(fscanf(fp , "%lf %lf¥n" , &t , &v)!=EOF)
	{
		Nol++;
	}
	rewind(fp);
	double *data = new double[Nol];
	for(int i=0;i<Nol;i++)
	{
		fscanf(fp , "%lf %lf¥n" , &t , &v);
		data[i] = v;
	}
	fclose(fp);

	int N = (int)(S/f+0.5);//入力信号の最短周期
	int M = Nol / N - 2;//M周期使ってfft

	fftw_complex *data1 , *data2;
	data1 = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*N*M);
	data2 = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*N*M);

	fp = fopen("data.txt" , "w");

	for(int i=N,n=0;n<N*M;i++,n++)//シミュレーション開始直後の値(~N)は捨てる
	{
		data1[n][0] = data[i];
		data1[n][1] = 0;
		fprintf(fp , "%e¥n" , data[i]);
	}
	fclose(fp);

	fftw_plan p = fftw_plan_dft_1d(N*M , data1 , data2 , FFTW_FORWARD , FFTW_ESTIMATE);
	fftw_execute(p);


	double sum = 0;
	for(int i=2*M;i<N*M/2+1;i+=M)//i=Mが入力信号と同じ周波数。i=M,2M,3M...を取り出す 
	{
		sum += data2[i][0]*data2[i][0]+data2[i][1]*data2[i][1];
	}
	printf("THD=%e¥n" , sqrt(sum) /
	sqrt(data2[M][0]*data2[M][0]+data2[M][1]*data2[M][1]));

 	delete[] data;
 	fftw_free(data1);
 	fftw_free(data2);

	return 0;
}
ジャジャガッチ | 電子工作 | 03:00 | comments(0) | trackbacks(0) |

スポンサーサイト

スポンサードリンク | - | 03:00 | - | - |
Comment









Trackback
URL:

07
--
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
--
>>
<<
--
PR
RECOMMEND
RECENT COMMENT
MOBILE
qrcode
OTHERS
Since 2013/09/17
LATEST ENTRY
CATEGORY
ARCHIVE
LINKS
PROFILE
SEARCH