• 締切済み

scilabのエラーについて

よろしくお願いします。 scilabにてオーバーラップを掛けてFFT解析をしたく、いくつかのWebを参考にさせて頂きながら 試しているのですが、”添字に誤りがあります”とエラーになってしまいます。 ソースは clear all; stacksize('max') v0=read_csv('C:\Users\*****\Desktop\scilab\data7.csv',",") v1 = evstr(v0);//上記v0データは文字列認識の為数値データに変換 //窓関数 N1=256; win_l = window('hn',N1); win_l_minus = win_l * -1; //FFTを実行 N=size(v1,1);//データ長 N2=N1/4; //N/2で50%オーバーラップ L = floor(N/N2)-1;//全サンプルに対するFFTの必要回数  y = 0;//出力ベクトルyの初期値 past_tail=zeros(1,N2);//ハーフオーバーラップ加算信号の初期値 for k=1:L n =N2*(k-1);//FFTの開始地点を更新   v2=win_l*v1( n+1 : n+N1 );//観測信号に窓を掛ける v3=fft(v2,-1);//FFT y1=fft(v3,1);//IFFT y=(y,past_tail + y1(1:N2));       //今回のIFFTの前半を前回のIFFT後半に加算してから出力ベクトルyに追加 past_tail=y1(N2+1:N1);//今回のIFFT結果の後半を記憶 end fv1dBuV=zeros(N,2); fv1dBuV(1,1) = 0;//DC =0Hzを入力 fv1dBuV(1,2) = 20*log10(abs(y(1)/N)/10.^-6);//dBuVに変換 sr = 4000; //sr:サンプリングレート[sample/sec] dt = 1/sr; // dt: サンプリング間隔[sec] T = N*dt; // T: 測定時間[sec] df =1/T; //Hz df: 1/測定時間->最低周波数、周波数分解能[Hz] for i=2:N fv1dBuV(i,1)=(i-1)*df///10^6;//周波数をMHzで記録 fv1dBuV(i,2)=20*log10(abs((y(i))/(N/2))/10.^-6); //N/2の左半分を使用して算出 end になります。(この後はPlot処理) y=(y,past_tail + y1(1:N2))で上記エラーが指摘されます。 プログラミングは初心者なのでどのように修正したらよいか 教えていただけると助かります。 また、オーバーラップを掛けてのFFT解析で他の方法があれば 教えてください。

みんなの回答

  • kmee
  • ベストアンサー率55% (1857/3366)
回答No.1

> ”添字に誤りがあります” と言われているのですから、添字が正しいかどうかを中心に見直すことです。 > y=(y,past_tail + y1(1:N2))で上記エラーが指摘されます。 ということなら、ここで使われている変数等に問題がある、ということです。 となると y pass_tail y1(1:N2) が関係していると考えます。 それぞれを見ると pass_tail=zeros(1,N2) y1=fft(v3,1);//IFFT です。 https://help.scilab.org/docs/5.5.0/ja_JP/zeros.html によると、 zeros(1,N2) は1×N2の行列です。 https://help.scilab.org/docs/5.5.0/ja_JP/fft.html によると、 fftが返すのはベクトルです。 「行列+ベクトル」って計算できましたっけ?

123nmnm
質問者

お礼

回答ありがとうございます。 上記エラー部分ですが、 http://www.kumikomi.net/archives/2010/09/ep30rir2.php?page=4 こちらを参考にさせていただいていたのですが。 今回ご指摘頂いた部分を見直してみます。 『行列+ベクトル』⇒修正してみます。

関連するQ&A

  • 複素数の絶対値について

    このプログラムでx2[i]y2[i]の絶対値の2乗がほしいんですけどどうすればよろしいで しょうか?? #include <stdio.h> #include <stdlib.h> #include <math.h> #define PI 3.14159265358979323846 /*省略*/ /* FFT main routine */ int fft(n, x, y) /*省略*/ /*** Main Program ***/ #define N 256 int main(argc, argv) int argc; char *argv[]; { FILE *fp; int i; static float x1[N], y1[N], x2[N], y2[N], x3[N], y3[N]; if ((fp = fopen(c1.txt, "r")) ==NULL){ printf ("File Open Error \"%s \" file.\n", c1.txt); exit(1); } for (i = 0; i < N; i++) { fscanf(fp, "%f", &x1[i]); x2[i] = x1[i] ; y1[i] = y2[i] = 0; } if (fft(N, x2, y2)) return EXIT_FAILURE; for (i = 0; i < N; i++) { x3[i] = x2[i]; y3[i] = y2[i]; } if (fft(-N, x3, y3)) return EXIT_FAILURE; printf(" Original data Fourier Transformed Inverse Transformed\n"); for (i = 0; i < N; i++) printf("%4d | %6.3f %6.3f | %6.3f %6.3f | %6.3f %6.3f\n", i, x1[i], y1[i], x2[i], y2[i], x3[i], y3[i]); return EXIT_SUCCESS; }

  • 配列のエラーに関して

    java言語を用いて,Householder変換を用いた固有値の数値計算に挑戦してみました.しかし,次のようなエラーが発生し上手くいきません.どなたかこの問題を解決するためにお力をかしていただけないでしょうか. ----------エラー内容-------------------------------------------------------------------------------- Exception in thread "main" java.lang.ArrayIndexOutOfBoundsException: Index 0 out of bounds for length 0 at Out.Mhouse(House.java:90) at House.main(House.java:10) ---------------------------------------------------------------------------------------------------- //Householder変換 public class House{ public static void main(String[] args){ double[][] A = new double[3][3]; int n = A.length; Out out = new Out(); for(int i = 0;i < n;i++){ for(int j = 0;j < n;j++){ if(j < n-1){ System.out.print(out.Mhouse(A)[i][j] + " "); }else if (j == n-1) System.out.println(out.Mhouse(A)[i][j]); }; }; }; }; class Out{ double[][] outpro(double[] x){ int n; n = x.length; double[][] A = new double[n][n]; for(int i = 0;i < n;i++ ){ for(int j = 0;j < n;j++){ A[i][j] = x[i] * x[j]; } } return A; }; double[][] Msca(double a,double[][] A){ int n = A.length; for(int i = 0;i < n; i++){ for(int j = 0;j < n;j++){ A[i][j] = a * A[i][j]; } } return A; }; double selfpro(double[] x){ double a = 0; int n = x.length; for(int i = 0;i < n; i++){ a = a + x[i] * x[i]; }; return a; }; double[] minus(double[] x, double[] y){ int n = x.length; double[] z = new double[n]; for(int i = 0;i < n;i++){ z[i] = x[i] - y[i]; }; return z; }; double[][] house_1(double[] x){ int n = x.length; double[][] A = new double[n][n]; for(int i=0;i < n;i++){ for(int j = 0;j < n;j++){ if(i == j){ A[i][j] = 1 - Msca(2/selfpro(x),outpro(x))[i][j]; }else{ A[i][j] = - Msca(2/selfpro(x),outpro(x))[i][j]; }; }; }; return A; }; double[][] house_2(double[] x){ double[][] z = new double[1][1]; z[1][1] = 1 - 2; return z; }; double[][] Mhouse(double[][] A){ int n = A.length; double[][] H = new double[n][n]; for(int i = 0;i < n;i++){ double[] x = new double[n-i]; double[] y = new double[n-i]; double[][][] L = new double[i][n-i][n-i]; for(int j = 0;j < n-i;j++){ x[j] = A[i][i+j]; if(j == 0){ y[j] = 1; }else{ y[j] = 0; }; x[j] = y[j] - x[j]; }; if(i < n-1){ L[i] = house_1(x); for(int k = 0;k < n-i;k++){ for(int l = 0;l < n-i;l++){ H[i+k][i+l] = L[i][k][l]; }; }; }else if(i == n-1){ L[i] = house_2(x); for(int k = 0;k < n-i;k++){ for(int l = 0;l < n-i;l++){ H[i+k][i+l] = L[i][k][l]; }; }; }; }; double[][] B = new double[n][n]; for(int i = 0;i < n;i++){ for(int j = 0;j < n;j++){ for(int k = 0;k < n;k++){ B[i][j] = H[i][k] * A[k][j]; }; }; }; return A; }; };

    • ベストアンサー
    • Java
  • Rubyのルンゲクッタ法がうまくいきません

    Rubyでルンゲクッタ法でy=exp(x**2)と比較しようという問題で、下のようなプログラムを組むとyの値がかなり大きくなってしまいます。どこが間違っているのでしょうか。教えてください。 #! ruby -Ks #ルンゲックッタ def df(a,b) z=2*a*b return z end n=10 x=Array.new(n) y=Array.new(n) x[0]=0.0 y[0]=1.0 h=1.0/n #任意で変更 fw=File.open("output2.txt","w") for i in 0..n k1=df(x[i],y[i])*h k2=df(x[i]+h/2,y[i]+k1/2)*h k3=df(x[i]+h/2,y[i]+k2/2)*h k4=df(x[i]+h,y[i]+k3) k=(k1+2*k2+2*k3+k4)/6 x[i+1]=x[i]+h y[i+1]=y[i]+k fw.print x[i]," ",y[i]," ","\n" end fw.close()

    • ベストアンサー
    • Ruby
  • Scilabを使ったジュリア集合の描画プログラム

    今、Scilabを使用してジュリア集合(充填およびそれ以外を含む集合)を描画するプログラムを書いています。 以前書いたC言語のプログラムをもとに書いているのですが、正確に描画できません。どうしたらよいでしょうか。教えてください。 与える条件は右上と左下の座標(複素数形式)と定数Cの値です。 以下に掲載したのが製作したプログラムです。 よろしくお願い致します。 //描画エリアの右上と左下の座標を複素数で設定する Z =[ -1.2-1.2*%i;1.2+1.2*%i]; //複素定数(C)を設定する C=0+0*%i; //描画エリアのx座標とy座標の各最小値と最大値を計算する。 xmin = min(real(Z)); xmax = max(real(Z)); ymin = min(imag(Z)); ymax = max(imag(Z)); Cr = real(C); Ci = imag(C); //描画点数を800×800に設定する。 N = 800; //各増分を計算する。 dx = (xmax-xmin)/(N-1); dy = (ymax-ymin)/(1-N); //プロットデータを"0"で初期化 map=zeros(N,N); //ジュリア集合の描画 i=1; for X=xmin:dx:xmax j=1; for Y=ymax:dy:ymin for k=1:30 x = X ^ 2 - Y ^ 2 + Cr; y = 2 * X * Y +Ci; if x^2 +y^2 > 4 then break; end map(j,i)=k; X=x; Y=y; end j=j+1; end i=i+1; end //プロットするための設定 Re = xmin:dx:xmax; Im = ymax:dy:ymin; clf(0); Sgrayplot(Re,Im,map');

  • エラーは出ませんが、実行結果ができません。

    このプログラムなんですが、エラーは出ませんが結果が 0群の項目1の正解率は0.000000です 1群の項目1の正解率は0.000000です 2群の項目1の正解率は0.000000です… この様になり、正解率がでません… 初心者で、わからないので困っています。 お願いします。 #include <stdio.h> #include <process.h> #define S 256 #define I 100 #define J 100 #define K 3 //グループの数 //#define M 50//サブコンテンツの数 void sum(int u[][J],int N,int n); void sort(int y[],int N,int u[][J],int n); void gunwake(int y[],int start,int N,int gunnum); void passege(int y[],int div[],int N,int num[],int u[][J],int n); static int y[I]; int div[K-1]; int divyouso=0; void main (void) { FILE *fp; int N=0,i=0,j=1,kou=0,n; //N:人数 n:問題数 static int u[I][J]; static int num[I]; char buf[S]; //ファイルオープン if ((fp=fopen("data_i2_1.csv","r"))==NULL){ printf("Can't open File\n"); exit(1); } // 問題数のカウント fgets(buf,S,fp); N+=1; while(buf[i]!='\n'){ kou=kou++; i+=1; } for(i=0;i<=kou;i=i+2){ u[N][j]=buf[i]-'0'; j=j++; } n=kou/2+1; // レコードの読み込み while (fgets(buf,256,fp)!=NULL){ N+=1; // 文字型から数値型へ変換 j=1; for(i=0;i<=kou;i=i+2){ u[N][j]=buf[i]-'0'; j=j++; } } sum(u,N,n); gunwake(y,0,N,K); passege(y,div,N,num,u,n); fclose(fp); } void sum(int u[][J],int N,int n) { //static int y[I]; int i,ii; //学習者iの得点の初期化 for(i=0;i<=I;i++) y[i]=0; //学習者iの得点の計算 for(i=1;i<=N;i++){ for(ii=1;ii<=n;ii++){ y[i]+=u[i][ii]; } } sort(y,N,u,n); } void sort(int y[],int N,int u[][J],int n) { int left,right,i,shift,t,v; static int num[I]; //学習者の番号記憶用変数numの初期化 for(i=0;i<=I;i++) num[i]=0; for(i=1;i<=N;i++) num[i]=i; //シェーカーソート left=0; right=N; while (left<right){ for(i=left;i<right;i++){ if(y[i]>y[i+1]){ t=y[i]; v=num[i]; y[i]=y[i+1]; num[i]=num[i+1]; y[i+1]=t; num[i+1]=v; shift=i; } } right=shift; for(i=right;i>left;i--){ if(y[i]<y[i-1]){ t=y[i]; v=num[i]; y[i]=y[i-1]; num[i]=num[i-1]; y[i-1]=t; num[i-1]=v; shift=i; } } left=shift; } } void gunwake(int y[],int start,int N,int gunnum){ int tmp; int i,up,down,real; if(gunnum>1){ tmp=N/gunnum+start; //printf("tmp:%d\n",tmp); for(i=tmp;y[tmp]==y[i];i--){ } down = i + 1; //printf("down:%d\n",down); for(i=tmp;y[tmp]==y[i];i++){ } up =i; //printf("up:%d/n",up); if(tmp-down > up-tmp) real=up; else real=down; div[divyouso]=real; divyouso++; printf("%d\n",real); gunwake(y,real,N-real,gunnum-1); } } void passege(int y[],int div[],int N,int num[],int u[][J],int n){ int div2[K+1]; int k=0,j,i; int pp[I][J]; div2[0]=0; div2[K]=n; for(i=0;i<K-1;i++){ div2[i+1]=div[i]; } for(k=0;k<K;k++){ for(j=0;j<n;j++){ pp[k][j]=0; for(i=div2[k];i<div2[k+1];i++){ pp[k][j]=pp[k][j]+u[num[i]][j]; } } } //確認 putchar('\n'); for(j=1;j<=n;j++){ for(k=0;k<K;k++){ printf("%d群の項目%dの正解率は%fです\n",k,j,pp[k][j]); } } }

  • Cプログラムによる画像の高速フーリエ変換FFT

    Cプログラムにる画像処理について教えてください。逆高速フーリエ変換(IFFT)によって周波数領域から実空間領域に変換する際、周波数帯域を(低周波数領域に)制限してIFFTを実行できるという文献があるのですが可能でしょうか?通常の逆離散フーリエ変換(IDFT)では可能でした。 以下のコードの修正で可能でしょうか? 例) 512*512の空間周波数領域データを212*212(仮に)の低周波数領域に帯域制限してIFFTで画像に戻したいのです。 以下の以下の一次元FFTコードを用いてX方向とY方向にIFFTしたいと考えています。 void FFT(int ir, int nx, float *xr, float *xi, float *si, float *co, unsigned short *brv) // 1次元フーリエ変換 // int ir; 順変換(1)と逆変換(-1) // int nx; 1次元FFTのデータ数 // float *xr; 実部のデータ xr[nx] // float *xi; 虚部のデータ xi[nx] // float *si; FFT用のサインデータ si[nx/2] // float *co; FFT用のコサインデータ co[nx/2] // unsigned short *brv; FFT用の入れ替えデータ brv[nx] { int i, j, n1, n2=nx, j3, j4, k, l, ll, d=1, g; float a, b, c, s; for(l = 1; l <= nx/2; l *= 2, d += d) { g = 0; ll = n2; n2 /= 2; for(k = 1; k <= n2; k++) { n1 = k-ll; c = co[g]; s = -ir*si[g]; g += d; for(j = ll; j <= nx; j += ll) { j3 = j+n1-1; j4 = j3+n2; a = xr[j3]-xr[j4]; b = xi[j3]-xi[j4]; xr[j3] += xr[j4]; xi[j3] += xi[j4]; xr[j4] = c*a+s*b; xi[j4] = c*b-s*a; } } } 通常のDFTでの帯域制限はこのように行い可能でした。 void fourier1d(int ir, float *fr, float *fi, int nx) { int i, j, n = 1; float *gr, *gi; double u, x; gr = (float *)malloc((unsigned long)nx*sizeof(float)); gi = (float *)malloc((unsigned long)nx*sizeof(float)); for(i = 0 ; i < nx ; i++) { u = i-nx/2; gr[i] = gi[i] = 0; for(j = 150 ; j < 362 ; j++) { x = j-nx/2; gr[i] += (float)( fr[j]*cos(2*PI*u*x/nx)+ir*fi[j]*sin(2*PI*u*x/nx)); gi[i] += (float)(-ir*fr[j]*sin(2*PI*u*x/nx)+fi[j]*cos(2*PI*u*x/nx)); } } if(ir == -1) n = 212; // 逆変換はデータ数で割る for(i = 0 ; i < nx ; i++) { fr[i] = gr[i]/n; fi[i] = gi[i]/n; } free(gr); free(gi); }

  • 2分探索法 『成功・不成功探索一回の実行時間を求める』

    C言語の授業で課題が出たのですが、自分が思っているような値が出なくて困っています。よろしければヒントだけでも頂ければなと思います。 ----------課題---------- 整列されているN個のデータに対して、その中にvがあるかどうかを判断する2分探索プログラムを実行する Nの値を10^6、5×10^6、10^7、5×10^7、10^8 と変化させたとき、成功探索、不成功探索一回の実行時間をそれぞれ求めよ。 このとき、Nと実行時間の関係はどのようになっているか(100字程度で)   N  成功探索  不成功探索  10^6  ***秒   ***秒 5×10^6  ***秒   ***秒  10^7  ***秒   ***秒 5×10^7  ***秒   ***秒  10^8  ***秒   ***秒 ---------------------------------------- ~↑を元に作成したプログラム(成功探索の場合)~ #include <stdio.h> #include <stdlib.h> #include <time.h> #define T 1000000 #define N 100000000 int a[N]; int main(void) { clock_t t1,t2; int t; int i; int l, r, k, v, m; for (i = 0; i < N; i++) a[i] = i; t1=clock(); for (t = 0; t < T; t++) { v = rand() % N; l=0; r=N-1; k=-1; while (r >= l) { m = (l+r)/2; if (v == a[m]) { k=m; break; } if (v < a[m]) r = m-1; else l = m+1; } } t2=clock(); printf("cpu time=%10.6f[micro sec]\n",(double)(t2-t1)*1000000/(CLOCKS_PER_SEC*T)); return 0; } ↑a[i]すべてに0~N-1を代入し、vにランダムに0~N-1の値を代入する。2分探索法で、v=a[i]となったら終了する。 実行時間は、↑の操作を10^6回行い、その平均を実行時間をする…単位:マイクロ秒 として、コンパイルして動いたんですが実行時間の値にずれがありどの値が一番適切か分からなくて困っています。 ↑のプログラムをそれぞれのNで実行したところ N=  10^6で実行・・・1.016 0.891 0.969 1.155 1.14 [micro sec] 5×10^6で実行・・・1.422 1.500 1.563 1.250 1.406 1.297  10^7で実行・・・1.750 1.360 1.37 1.407 1.672 1.531 5×10^7で実行・・・1.859 1.797 1.812 1.800  10^8で実行・・・2.062 2.140 2.320 2.500 このような結果になりました。 これで正しいのでしょうか?もう少しずれの幅が小さいと決められそうなのですが…そもそもプログラムが間違ってるんでしょうか? でも、Nが大きくなるにつれて実行時間が増えてるのでこちらはまだいいんですが・・・問題は不成功探索の方です。 次に 不成功探索では↑のプログラムの 乱数vのところを変化させて v=rand() % N;という箇所を v=N;としました。 a[i]には0~N-1が入っているので、v=Nとすれば必ず不成功になるはずですよね? こうして実行してみると N= 10^6、5×10^6、10^7、5×10^7、10^8と値を変化させても 0.719 0.54 0.625 0.534 [micro sec] に近い値ばかり出てしまい、正しい値とは到底思えません。 不成功探索は成功探索より時間がかかるはずですよね?なのになぜこのようになってしまったのでしょうか? 後、Nと実行時間の関係とは、最終的に得られた結果を元にして「実行時間はlog2Nに比例している」と言えばいいんでしょうか? こういうものってどのように回答したらいいのかヒントだけでももらえると非常に有難いです。 長々とスイマセン。 少し気づいたことなど些細なことでも全然かまわないので、どうかよろしくお願いします。

  • 物理の問題です。よろしくお願いします。

    物理の問題です。よろしくお願いします。 環状鉄心の左右にコイルN1(60回巻き)、N2(30回巻き)が巻かれている。 N1について、0.1[sec]間に3[A]の割合で変化させた場合、 コイルN1にはeが150[V]誘起された。 このとき、相互インダクタンスMを求めよ。 という問題ですが、 自己インダクタンスL1はe=L1(ΔI/Δt)より、5[H]とまでは分かるのですが、 M=L1(N2/N1)=5(30/60)=2.5[H]となる過程が理解できません。 ご教授願います。

  • 関数の出力引数をクラスにするには?

    既出、または基礎の質問でしたらすみません。 ここでも他の検索エンジンでも見つけられなかったので。。。 C++です。 クラスを出力する関数を作りたいのですが、うまくできません。 ソースは以下のとおりです。問題は、プログラム下方のf1(),f2(),main()です。 長くて、そして見づらくてすみません・・・ //////////// #include<stdio.h> class test{ private: int num; float *vec; public: test(int n=1); //ctor ~test(); //dtor int getnum(){return num;} float* getvec(); void set(int,float*); void show(); }; test::test(int n){ num = n; vec = new float[n]; for(int i=0; i<n; i++) vec[i] = (float)i; } test::~test() {delete[] vec;} float* test::getvec(){ float *v; v = new float[num]; for(int i=0; i<num; i++) v[i] = vec[i]; return v; } void test::set(int n, float *v){ num = n; vec = new float[n]; for(int i=0; i<n; i++) vec[i] = v[i]; } void test::show(){ for(int i=0; i<num; i++) printf("%d: %g\n",i,vec[i]); } void f1(test &x, test &y){ int n; float *v; n = x.getnum(); v = x.getvec(); for(int i=0; i<n-1; i++) v[i] = 2.0*v[i]; y.set(n-1,v); } test f2(test x){ test y; //* int n; float *v; n = x.getnum(); v = x.getvec(); for(int i=0; i<n-1; i++) v[i] = 2.0*v[i]; y.set(n-1,v); return y; //** } void main(){ test x,y; int n = 4; float v[4] = {1.0,2.0,3.0,4.0}; printf("x:\n"); x.set(n,v); x.show(); printf("f1:\n"); f1(x,y); y.show(); printf("f2:\n"); y=f2(x); y.show(); } //////////// これを実行すると x: 0: 1 1: 2 2: 3 3: 4 f1: 0: 2 1: 4 2: 6 f2: 0: 7.38979e-38 1: 7.38979e-38 2: 6 となります。 関数f2がうまく動かない理由がわかりません。。。 出力引数にクラスはとれないのでしょうか?? よろしくお願いします。

  • 行列と数列の関係式に関する問題(立教)

    これも今年の立教大学の問題です。特に(iii)以降、よろしくお願いします。   l 7 18 | A= | -3 -8 | とおく.Aに対して        | x[1] |  | x[n+1] |   | x[n] |        | y[1] |, | x[n+1] | = A| y[n] | により座標平面上の点P[n](x[n],y[n])(n=1,2,…)を定める. このとき, 次の問(i)~ (iv)に答えよ. (i) P[2], P[3] の座標を求めよ. (ii) すべての自然数 n について, P[n]が座標平面上のあるひとつの直線 l 上にあるこ  とを示せ また, 直線 l の方程式を求めよ. (iii) x[n+1] を x[n] の式で表せ. (iv) x[n],y[n] を n の式で表せ