• 締切済み

カーブフィットを行うプログラムについて

カーブフィットを行うプログラムについて質問があります。 下記は、周波数応答関数を入力し、 カーブフィットを行うプログラムですが、 プログラムの下4行の意味が全く分かりません。 (下4行とは、■マークが先頭に書いてあるものです) どうか何卒ご教授お願い致します。 わかりにくい質問の仕方をしてしまい、申し訳ないです。 ---------------------- clear variables global;comgui('close all') fm='wagawara_ryou10_'; n1=101; n2=142; n=0; for i=n1:n2 fnm=strcat(fm,int2str(i),'_h.txt') [IIw,Mag,Phi,Q,Q]=textread(fnm,'%f %f %f %f %f','headerlines',3); n=n+1; IIxf(:,n)=complex(Mag.*cos(Phi*pi/180),Mag.*sin(Phi*pi/180)); end %---初期固有振動数データと減衰比データの読み込み.2005.9.13--- fn='wagawara_ryou10_mif_pk_dm.txt' [peak,zeta]=textread(fn,'%f %f','headerlines',2); pkze=[peak,zeta]; %---処理開始--- iiplot %伝達関数のプロット XF(5).po=pkze %初期固有振動数と減衰比を代入。 idcom('e 15 371'); %idcom('e i w') i:帯域幅,w:中心周波数 idcom('est'); idcom('eup .05 .002 -10'); %idcom('eup dstep fstep num i') ■out=id_rm(XF(5),[1 1 1 1]); ■XF(3)=res2xf(out,IIw); ■iicom('IIxhOn'); ■[som2,ga2,pbs2,cps2]=res2nor(IIres,IIpo,IDopt); ---------------------- よろしくお願い致します。 失礼します。

みんなの回答

  • noocyte
  • ベストアンサー率58% (171/291)
回答No.4

「+MATLAB +"Structural Dynamics Toolbox"」で Google 検索 http://www.google.co.jp/search?sourceid=navclient-ff&ie=UTF-8&rls=GGGL,GGGL:2006-34,GGGL:ja&q=%2BMATLAB+%2B%22Structural+Dynamics+Toolbox%22 してみましたが,日本語ではほとんど情報がないですね. 次のものが多少は参考になるかも. ANSYSとMATLABを用いた構造系と制御系の同時最適化 ‐スマート構造への適用‐ (2006/03/03 日本機械学会「設計研究会」発表スライド) http://www.jsme.or.jp/dsd/A-TS12-05/minutes/18/Ishizuka_2006-03-03ANSYS_MATLAB.pdf Structural Dynamics Toolbox 5 等のカタログ http://www.sdtools.com/pdf/datasheet_japan.pdf MATLAB:適用分野:Test & Measurement:アプリケーション例 (サイバネットシステム) Structural Dynamics Toolboxを使用した振動解析/制御 http://www.cybernet.co.jp/matlab/solution/applications/t_m/application1.shtml

yasuyasu19
質問者

お礼

noocyteさんへ 本当に、色々な情報をありがとうございます。 凄くありがたいです。 教えていただいた情報を無駄にしないよう、 一生懸命勉強いたします。 ありがとうございました。 失礼します。

  • noocyte
  • ベストアンサー率58% (171/291)
回答No.3

#1 です. 「+MATLAB +id_rm」で Google 検索 http://www.google.co.jp/search?sourceid=navclient-ff&ie=UTF-8&rls=GGGL,GGGL:2006-34,GGGL:ja&q=%2BMATLAB+%2Bid_rm してみると,id_rm のマニュアル (後述) と,#2 さんのおっしゃっていた Structural Dynamics Toolbox のマニュアルが見つかりました. Structural Dynamics Toolbox FEMLink (For Use with MATLAB(R)) http://aertia.com/docs/sdtools/sdt_manual.pdf んで,色々検索しまくって,■の行にある id_rm のマニュアルの一部だけ訳してみました. (合ってるかどうかわかりません.制御工学なんて20数年ぶりですし,こんなに検索しまくったのは初めてです.) ●OUT = id_rm(IN, multi) http://www.sdtools.com/help/id_rm.html 用途:MIMO (多入力多出力) 制御系の最小モデルを作成し,    scaled modal 入出力を得るために reciprocity (相互?) 制約を適用する. 入力:(1) IN:必要なフィールドは下記.      IN.res:留数.      IN.po:極.      IN.idopt:同定オプション.       .FittingModel:Posit,Complex,Normal モード.       .NSNA:センサ/アクチュエータの数.       .Reciprocity:未使用,1 FRF (周波数応答関数?) または真の MIMO.       .Collocated:Reciprocity を使用する場合,collocated FRF の        indices (指数?).    (2) multi:IN.po 内の各極の (重根の) 多重度を表すベクトル (オプション). 出力:OUT:下記のフィールドを持つ.(将来変更されそう.)     .po:(適切な多重度を持つ) 極.     .def:出力 shape (シェープ?) 行列 (CPSI).     .DOF:センサの DOF (自由度?)     .psib:入力 shape 行列 (PSIB).     .CDOF:collocated FRF の指数.     .header:ヘッダ (最大72文字のテキスト×5行). その他の組み込み関数は, ・res2xf:留数モデルに関連付けられた多項式表現を作成する.  留数と極の集合に対し,対応する FRF (周波数ポイントwで評価される) を生成する. ・res2nor:複素留数を,ノーマルモード留数または比例減衰ノーマルモード形式に近似変換する. ・XF:Database wrapper object のためのユーザ・インタフェース. ・iicom:FRF データ可視化のための UI (ユーザ・インタフェース) コマンド関数.  yasuyasu19 さんが,振動モード解析 (または制御工学) と MATLAB について どの程度の知識をお持ちなのかわかりませんが,前者について充分な知識がおありなら, 上記の Structural Dynamics Toolbox のマニュアルを読めばわかると思います. その場合,MATLAB についてわからないことは検索するなり,ここで質問すれば いいと思います.  しかしそうでない場合 (例えば上の訳を見てもピンとこない,「留数」や「極の多重度」 の意味がわからない) には,ここで誰かに (私は無理です orz) 教えてもらったり,検索 したりしたとしても,失礼ながら (この分野を勉強しない限り) どだい無理です.職場に この専門分野に詳しい人がいるはずですので,その人に手取り足取り教えてもらうか, さもなくば「制御工学がわからないのでできません」宣言をして他の人に代わってもらう しかないと思います. ●参考 (になるかもしれない) 「"モーダルパラメータ" "カーブフィット"」で Google 検索 http://www.google.co.jp/search?sourceid=navclient-ff&ie=UTF-8&rls=GGGL,GGGL:2006-34,GGGL:ja&q=%22%E3%83%A2%E3%83%BC%E3%83%80%E3%83%AB%E3%83%91%E3%83%A9%E3%83%A1%E3%83%BC%E3%82%BF%22+%22%E3%82%AB%E3%83%BC%E3%83%96%E3%83%95%E3%82%A3%E3%83%83%E3%83%88%22 モーダル解析技術 (教材情報) http://www.tetras.uitec.ehdo.go.jp/db/kyouzai/kyouzai_detail.php?stat=3&id=177 振動モード解析の理論と実際 (教材情報) http://www.tetras.uitec.ehdo.go.jp/db/kyouzai/kyouzai_detail.php?stat=2&id=52 古典制御へのMATLABの利用 http://ichiro.maruta.googlepages.com/matlabtutorial

yasuyasu19
質問者

お礼

こんにちは! noocyteさん、色々と調べて頂き、本当に感謝します。 教えて頂いたHPをもとに、調べてみようと思います。 再度質問をするかもしれません。 その際は、どうぞよろしくお願い致します。 本当にありがとうございました。 失礼します。

  • rabbit_cat
  • ベストアンサー率40% (829/2062)
回答No.2

matlabのスクリプトでしょうか。 見慣れない関数は、Structural Dyanamics Toolboxのもののようですが、私は使ったことがないのでよくわかりません。

  • noocyte
  • ベストアンサー率58% (171/291)
回答No.1

> どうか何卒ご教授お願い致します。 それは無理です.理由は (少なくとも) 2つ. (1) この言語が何という言語なのか示されていません.   だからこの言語を知らない人には構文の意味さえわからないし,   調べようもありません.   (たまたまこの言語を使っている人がここを見ていれば,    何かコメントがつくかもしれませんが.) (2) (特に■の行の) どれがユーザ定義関数で,どれが組み込み関数   なのかわかりません.ユーザ定義関数ならば,その定義も示して   いただかないと,それが何をしているかなんて誰にもわかりません.   組み込み関数ならば,(1) の理由で調べようがありません.

yasuyasu19
質問者

お礼

説明不足でした。 折角見ていただいたのに、本当に申し訳ございません。 noocyteさんの問いに対する回答を書きましたので どうかご覧下さい。 ---------- (1) この言語が何という言語なのか示されていません.   だからこの言語を知らない人には構文の意味さえわからないし,   調べようもありません. ↑↑ 言語はMATLABです。 (2) (特に■の行の) どれがユーザ定義関数で,どれが組み込み関数   なのかわかりません.ユーザ定義関数ならば,その定義も示して   いただかないと,それが何をしているかなんて誰にもわかりません. ↑↑ 組み込み関数と入力データを以下に示します。それ以外は、組み込み関数です。 'wagawara_ryou10_'←入力データ strcat←組み込み関数 textread←組み込み関数 complex←組み込み関数 'wagawara_ryou10_mif_pk_dm.txt'←入力データ textread←組み込み関数 idcom←組み込み関数 ■id_rm←組み込み関数 ■res2xf←組み込み関数 ■iicom←組み込み関数 ■res2nor←組み込み関数 ---------- 以上です。 返信が遅くなってしまい、大変申し訳ございませんでした。 まだまだ説明不足の箇所があるかと思いますが、 これに回答していただけると本当にありがたいです。 どうかよろしくお願い致します。 失礼します。

関連するQ&A

  • MATLABでカーブフィットしたデータをテキスト形式で保存する方法を教えてください

    下記は、周波数応答関数を入力し、 カーブフィットを行うプログラムです。 このプログラムを実行し、カーブフィットをできるのは いいのですが、どうすればそのデータをテキスト形式で保存し, エクセルなどで実行できるのかがわかりません。 どうか何卒ご教授お願い致します。 ---------------------- clear variables global;comgui('close all') fm='wagawara_ryou10_'; n1=101; n2=142; n=0; for i=n1:n2 fnm=strcat(fm,int2str(i),'_h.txt') [IIw,Mag,Phi,Q,Q]=textread(fnm,'%f %f %f %f %f','headerlines',3); n=n+1; IIxf(:,n)=complex(Mag.*cos(Phi*pi/180),Mag.*sin(Phi*pi/180)); end %---初期固有振動数データと減衰比データの読み込み.2005.9.13--- fn='wagawara_ryou10_mif_pk_dm.txt' [peak,zeta]=textread(fn,'%f %f','headerlines',2); pkze=[peak,zeta]; %---処理開始--- iiplot %伝達関数のプロット XF(5).po=pkze %初期固有振動数と減衰比を代入。 idcom('e 15 371'); %idcom('e i w') i:帯域幅,w:中心周波数 idcom('est'); idcom('eup .05 .002 -10'); %idcom('eup dstep fstep num i') out=id_rm(XF(5),[1 1 1 1]); XF(3)=res2xf(out,IIw); iicom('IIxhOn'); [som2,ga2,pbs2,cps2]=res2nor(IIres,IIpo,IDopt); ※注意 'wagawara_ryou10_'←入力データ strcat←組み込み関数 textread←組み込み関数 complex←組み込み関数 'wagawara_ryou10_mif_pk_dm.txt'←入力データ textread←組み込み関数 idcom←組み込み関数 id_rm←組み込み関数 res2xf←組み込み関数 iicom←組み込み関数 res2nor←組み込み関数 ---------------------- どうぞよろしくお願いします。

  • MATLABでカーブフィットしたデータを、関数「xlswrite」「save」を用いて、xls形式・txt形式で保存する方法を教えてください

    下記は、周波数応答関数を入力し、 カーブフィットを行うプログラムです。 カーブフィットのデータをエクセル形式、テキスト形式で 保存をする場合、関数[save]または[xlswrite]を使用し, どのように書けばいいのでしょうか。 どうか何卒ご教授お願い致します。 ---------------------- clear variables global;comgui('close all') fm='wagawara_ryou10_'; n1=101; n2=142; n=0; for i=n1:n2 fnm=strcat(fm,int2str(i),'_h.txt') [IIw,Mag,Phi,Q,Q]=textread(fnm,'%f %f %f %f %f','headerlines',3); n=n+1; IIxf(:,n)=complex(Mag.*cos(Phi*pi/180),Mag.*sin(Phi*pi/180)); end %---初期固有振動数データと減衰比データの読み込み.2005.9.13--- fn='wagawara_ryou10_mif_pk_dm.txt' [peak,zeta]=textread(fn,'%f %f','headerlines',2); pkze=[peak,zeta]; %---処理開始--- iiplot %伝達関数のプロット XF(5).po=pkze %初期固有振動数と減衰比を代入。 idcom('e 15 371'); %idcom('e i w') i:帯域幅,w:中心周波数 idcom('est'); idcom('eup .05 .002 -10'); %idcom('eup dstep fstep num i') out=id_rm(XF(5),[1 1 1 1]); XF(3)=res2xf(out,IIw); iicom('IIxhOn'); [som2,ga2,pbs2,cps2]=res2nor(IIres,IIpo,IDopt); ※注意 'wagawara_ryou10_'←入力データ strcat←組み込み関数 textread←組み込み関数 complex←組み込み関数 'wagawara_ryou10_mif_pk_dm.txt'←入力データ textread←組み込み関数 idcom←組み込み関数 id_rm←組み込み関数 res2xf←組み込み関数 iicom←組み込み関数 res2nor←組み込み関数 ---------------------- どうぞよろしくお願いします。

  • 周波数応答関数の計算

    打撃実験を行って解析するのですが下記のプログラムが分かりません。 [IIw,Mag,Phi,Q,Q]=textread(fnm,'%f %f %f %f %f','headerlines',3); IIxf(:,n)=complex(Mag.*cos(Phi*pi/180),Mag.*sin(Phi*pi/180)); 読み込んだファイルの内容が Np Nc / Freq Hm Ha Coh XPwr 6401 5 0 11.8941 0 0.938239 1.78076E-006 1 27.7519 130.402 0.915975 2.26936E-008 2 3.95322 -38.8943 0.973454 1.64854E-008 . . 一列目が周波数と思います。2列と3列が何か分かりません。 また、IIxfは周波数応答関数だと思うのですがどの公式を用いているのでしょうか?

  • c言語でDFTのプログラムを作成したのですが

    c言語でDFTのプログラムを作成しました。 以下にソースを載せます。 #include<stdio.h> #include<math.h> #include<stdlib.h> #include<time.h> #define PI 3.141592653589793 #define N 64 //データ数 DFT(double result[]){ int i,k; double A[N],B[N],T=0; //A[N]:実数部,B[N]:虚数部 double a,b; for(k=0;k<N;++k){ a=b=0; for(i=0;i<N;++i){ a=a+result[i]*cos(2*PI*i*k/N);      b=b+(-1.0)*result[i]*sin(2*PI*i*k/N); } A[k]=a/N; B[k]=b/N; } for(i=0;i<N;++i){ printf("[%f秒]:Re:%f,Im:%f\n",T,A[i],B[i]); //変換後の値を表示 T=T+(0.1/N); } } main(){ int i; double T=0; double Sampdata; double result[N]; for(i=0;i<N;++i){ Sampdata=5*sin(20*PI*T);      //0~0.1秒間をN個にサンプリング result[i]=Sampdata; //サンプリングデータを代入 T=T+(0.1/N); } clock_t start,end; //処理時間計測開始 start=clock(); DFT(result); end=clock(); printf("%.2f秒かかりました\n",(double)(end-start)/CLOCKS_PER_SEC); //処理時間表示 } 元信号には5sin(20πt)の値を入れています。この信号は周期は0.1secです。 これでフーリエ変換を行うとデータ数N/2を中心に対称なデータが出てくるのですが、処理が終わるのが早い気がするんです。 例えば2^15個のデータで実行しても2分もかからずに処理が終わってしまいます。一応、対称性が出てるとはいえ、終わるのが早すぎる気がするのですが、おかしい所があれば教えていただけると嬉しいです。 よろしくお願いします。

  • 漸化式のプログラム

    n=100としてn+1個の点(Xj,Yj) (j=1,2,3・・n)はどのようなグラフになるか?ただし、h=1/n, Xj=jh, Y(j+1)=Yj+hYj, Y(0)=1である。 このプログラムを教えてください。 #include <stdio.h> #include <math.h> main() { double pi, x, y; int i; pi=4*atan(1); x=1/100; y=1; for(i = 1; i<=100; ++i){ x=i*1/100; y[n+1]=y[n]+1/100*y[n]; printf("%lf %lf\n", i, x, y,); } } ではだめなんでしょうか?

  • C++言語のプログラムをfortranに変換

    const int N = 100; const double q = 10.0, dt = 0.00001, Dm = 30.0, t0 = 2.0, K = 1.0, pi = 3.141592, f = 3.0; double C[N], dC[N]; double dx = q/N; for (int i = 0; i < N; ++i) C[i] = 0; // 初期条件 for (double t = 0; t < t0; t += dt) {  C[0] = 0; // 境界条件1  C[N - 1] = K*sin(2*pi*f*t); // 境界条件2  for (int i = 1; i < N - 1; ++i) dC[i] = (Dm*(C[i + 1] - 2*C[i] + C[i - 1])/(dx*dx))*dt;  for (int i = 1; i < N - 1; ++i) C[i] += dC[i]; } for (int i = 0; i < N; ++i) cout << C[i] << endl; // t = t0 このプログラムをfortranに変換できる方いますか?

  • 離散フーリエ変換(DFT)のプログラム

    タイトルの通りのレポートを出されたのですが、その問題に似たサンプルソースすら何をいっているのかわからない状態です。ひとまず、サンプルソースが何をいっているのか理解したいので、いくつか教えてください。 ソースです。質問はその後に書かせてもらいました。 #include<stdio.h> #include<math.h> #define N 10 #define F 0.1 #define PI 3.14151692 #define SQR(x) ((x)*(x)) void func() { int n; FILE *fp; fp=fopen("temporal.data","w"); for(n=0;n<N;n++) fprintf(fp,"%lf\n",cos(2.0*PI*F*(double)n)); fclose(fp); } void get_data(double x[]) { int n; FILE *fp; fp=fopen("temporal.data","r"); for(n=0;n<N;n++) fscanf(fp,"%lf",&x[n]); fclose(fp); } void dft(double x[],double X_r[],double X_i[]) { int k,n; for(k=0;k<N;k++){ X_r[k]=X_i[k]=0.0; for(n=0;n<N;n++){ X_r[k]+=x[n]*cos(2.0*PI*(double)n*(double)k/(double)N); X_i[k]-=x[n]*sin(2.0*PI*(double)n*(double)k/(double)N); } } for(k=0;k<N;k++) printf("X[%d]=%lf+j%lf\n",k,X_r[k],X_i[k]); } void amplitude(double X_r[],double X_i[]) { int k; FILE *fp; double amp; fp=fopen("amp.data","w"); for(k=0;k<N;k++){ amp=sqrt(SQR(X_r[k])+SQR(X_i[k])); fprintf(fp,"%lf\n",amp); } fclose(fp); } main() { double x[N],X_r[N],X_i[N]; func(); get_data(x); dft(x,X_r,X_i); amplitude(X_r,X_i); } 文字数の制限があるみたいなので、質問を別にさせてもらいます。

  • ガウス関数を少しずつずらして足し上げるプログラム

    ある値をある検出器に入力したときに出力結果がその入力したを中心にガウス関数のように広がってしまう場合を考え,ある関数の値を小刻みに入力したとき、その関数が出力されたときそれぞれ値のガウス関数の広がりにより重ね合わされてどのように見えるか数値計算してみようと思い、以下のソースコードを書いて見ました。行列で考えてそれぞれの値の寄与を足しあげて見ようとしたのしたのですが、実行結果をみると明らかにおかしい10^19などの数値が見られます。これはなぜでしょうか。私のソースコードのどこに問題があるのでしょうか。お手数をお掛けしますが回答よろしくおねがいします。 インデントが反映されていませんでしたらすみません。 ---------------------------ソースコード------------------------------------------ #include<stdio.h> #include <stdlib.h> #include<math.h> #define N 100 #define f(x) pow(x,-1.7) //この関数が検出器に入力したときどのように見えるか知りたい// void mat_product(double ec[N],double r[N][N],double e[N]) { int i,j; for(i=0; i<N; i++){ for(j=0; j<N;j++){ ec[i]+=r[i][j]*e[j]; } } } int main(void) { int i,j; double ec[N],r[N][N],e[N]; double max,min, ei,E,sigma,min1,min2,m,pi; sigma=0.08; //keV// max=10.5; min=0.5; min1=0.5; min2=0.5; m=0.5; ei=(max-min)/N; pi=3.141592; for(i=0; i<N ;i++){ for(j=0; j<N ;j++){ r[j][i]=(1/(pow(2*pi,0.5)*sigma))*exp(-1*pow(m-min1,2)/(2*pow(sigma,2))); m+=ei; } e[i]=f(min2); m=0.5; min1+=ei; min2+=ei; } mat_product(ec,r,e); for(i=0; i<N; i++){ printf("%g\t%g\n",min,ec[i]); min+=ei; } return 0; } ---------------------------出力結果--------------------------------- 0.5 22.0512 0.6 23.8139 0.7 18.9 0.8 14.9642 0.9 12.1864 1 10.1512 1.1 8.60991 1.2 7.4112 1.3 6.45833 1.4 5.68688 1.5 5.05253 1.6 4.52387 1.7 4.07814 1.8 3.69846 1.9 3.37209 2 3.08926 2.1 2.84239 2.2 2.62547 2.3 2.43375 2.4 2.26337 2.5 2.1112 2.6 1.66592 2.7 1.85166 2.8 1.7404 2.9 1.63941 3 1.54742 3.1 1.46337 3.2 1.38635 3.3 1.31558 3.4 4.47979e+30 3.5 1.19017 3.6 1.13444 3.7 1.08275 3.8 1.0347 3.9 0.989949 4 0.9482 4.1 0.909181 4.2 0.872652 4.3 0.838401 4.4 0.806238 4.5 0.775991 4.6 0.747509 4.7 0.720652 4.8 0.695296 4.9 0.671329 5 0.648648 5.1 0.627161 5.2 0.606783 5.3 0.587437 5.4 0.569053 5.5 0.551566 5.6 0.534917 5.7 0.519053 5.8 0.503924 5.9 0.489483 6 0.475689 6.1 0.462502 6.2 0.449887 6.3 0.437809 6.4 0.426239 6.5 0.415146 6.6 0.404506 6.7 0.394292 6.8 0.384482 6.9 0.375054 7 0.365988 7.1 0.357265 7.2 0.348868 7.3 0.34078 7.4 0.332986 7.5 0.325471 7.6 0.318222 7.7 0.311226 7.8 0.304472 7.9 0.297947 8 0.291642 8.1 0.285546 8.2 0.27965 8.3 0.273945 8.4 -0.0404521 8.5 0.263075 8.6 0.257894 8.7 0.252874 8.8 0.248007 8.9 0.243288 9 0.238709 9.1 7.46629e+19 9.2 0.229953 9.3 0.225764 9.4 0.221696 9.5 0.217743 9.6 0.2139 9.7 0.210164 9.8 0.206531 9.9 0.202996 10 0.199557 10.1 0.196209 10.2 0.192868 10.3 0.185672 10.4 0.140717

  • 桁落ちのプログラムで真の値と計算結果

    #include <stdio.h> const double PI=3.141592653589793; double sum(long m) { double n,term,sum; n=1; sum=0; term= 1.0/ (n*n);/*初項*/ while( n<=m ){ sum+=term; n++; term= 1.0/ (n*n);/*次項の計算*/ } return sum; } /*この計算の答えはπ*/ int main(void) { double s; long m; int i; m=1; for(i=0;i<9;i++){ s=sum(m); printf("%2d m=%10ld sum= %22.16e err= %22.16e \n",i,m,s,s-PI*PI/6); m*=10;/*次は10倍にする*/ } return 0; で真の値と計算結果を調べるにはどうしたらいいのでしょうか?

  • 部分群の位数の形は?

    自然数n>1に於いて  m n=Π(pi)^(ri) (m,ri∈N, p1,p2,…,pmは相異なる素数)  i=1 とする。 位数nの群の部分群の位数は (pi)^(ri) (m,ri∈N,piは素数 (i=1,2,…,m)) のみである。 (つまり、(pi)^(ri)pj^(rj)といった形の位数を持つ事はない) という主張なのですがこの命題は正しいといっていいでしょうか?