• 締切済み

ゆらぎの分析

MATLABを使ってゆらぎの分析をしたいのですが、 FFT後、両対数軸で表示しようとすると、「警告:負のデータは無視されます」と出てきて、グラフが表示されなかったり、1本の直線が表示させたり、上手く表示できません。 FFTをする、対数軸で表示、の間にワンクッショなにかが必要なのでしょうか? 詳しい方教えてください。 一応テキストを見ながら作ったのですが、貼っておきます。 clear; fs=10000; dft_size=1024; x=wavread(' .wav');  //空けてあります。 w=HanningWindow_(dft_size); for n=1:dft_size, x(n)=x(n)*w(n); end X=fft(x,dft_size); A=zeros(1,dft_size/2+1); frequency=zeros(1,dft_size/2+1); for k=1:dft_size/2+1, A(k)=20*log10(abs(X(k))); frequency(k)=(k-1)*fs/dft_size; end subplot(2,1,1), plot(frequency,A); xlabel=('周波数[Hz]'); ylabel=('振幅[db]'); % 対数軸で表示 subplot(2,1,2), loglog(frequency,A); xlabel=('周波数[Hz]'); ylabel=('振幅[db]');

みんなの回答

回答No.1

単純にlog(-1)はxy平面には表現できないだけです。 log(0)=-∞もグラフには当然表示できません。 0以下だった場合0.001に置き換えるなどの処理があれば、 表示できます。

関連するQ&A

  • MATLABで二次元フーリエ変換

    画像処理のプログラムを作成しています。参考にしている参考書は【最新MATLABハンドブック】という本です。この本を参考にして、一次元フーリエ変換のプログラムから画像処理の二次元フーリエ変換のプログラムに変更させたいのですが、fft2のところでエラーが出てしまいます。自分なりにプログラムを書き直してみたのですが、fft→fft2に関数変更する以外にもっと根本的なことが必要なのでしょうか? 作成したプログラムはこれです。 clear;close all;n=256;dt=0.005; t=((1:n)-1)*dt; f=t/dt/dt/n;n2=n/2;n2p1=n2+1; X=imread('001.bmp'); X=rgb2gray(X); X=double(X); [m,n]=size(X); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %calc corresponding index number F=50; index=round(F*dt*n+1);index1=(index-1):(index+1); index2=n+2-index1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %band elminate filter befil=ones(size(X)); befil(index1)=zeros(size(index1)); befil(index2)=zeros(size(index2)); Y=fft2(X);    ←ここでエラーが出ます。 subplot(2,2,1) plot(f,abs(Y)/n2,'r');axis([0 f(n2p1) 0 1]) xlabel('Frequency (Hz)'),ylabel('Magnitude'); title('Original Signal in freq domain'); %%Now apply fft filter in freq domain fftX=fftX.*befil; subplot(2,2,3) plot(f,abs(Y)/n2,'r');axis([0 f(n2p1) 0 1]) xlabel('Frequency (Hz)'),ylabel('Magnitude'); title('filtered Signal in freq domain'); subplot(2,2,2) imaagesc(X),colormap(gray),axis tight; subplot(2,2,4) imagesc(real(ifft2(Y))),colormap(gray),axis tight; 画像の一部分の情報のみを欠落させたいのです。そのためにはバンドエルミネーションフィルタと思ったのですが…もし違うようならご指摘をお願いします。ちなみに001.bmpはカラー画像です。

  • パワースペクトルの0dBの値について

    matlabを用いて時系列データのパワースペクトルを計算、対数表示でグラフにプロットしています。 matlabのfft関数のヘルプに記載されていた例を参考に、以下のような自作の関数を組んで計算しているのですが、2つ質問があります。 1. この場合、パワースペクトルの0dBは、時系列の値ではいくつになるのでしょうか? 0dB = 1でしょうか? 2. 関数の中において、power_fftをN_fft(高速フーリエ変換時のデータ長さ)で割っている事の理由がわからなかったのですが、なぜ行っているのでしょうか? ------------------------------------------------ 以下、自作のパワースペクトル表示関数を示します function Power_db(Data_xx) fs = 1000; % サンプリング周波数 [Hz] N_fft = 2^12; % 4096個 xx_fft = fft(Data_xx,N_fft); power_xx = xx_fft.* conj(xx_fft) / N_fft; frequency = fs*(0:2047)/N_fft; % (0:2047)で500Hzまで plot(frequency,20*log10(power_xx(1:2048))) xlabel('Frequency [Hz]'); ylabel('Power Spectrum [dB]'); grid on; --------------------------------------------------------- 以上です。 何か足らない情報などがありましたらお申し付けください。 よろしくお願い致します。

  • MATLABのグラフで軸目盛りのフォントサイズを大きくする方法?

    MATLABのplotグラフで、 xlabel、ylabel、titleなどは、 FontSizeプロパティの指定により、 サイズを簡単に変更できるのですが、 軸の目盛りのフォントサイズを変更する方法が分からず 困っています。どなたか、キーワード、参考URLなど ご教示頂ければありがたいです。よろしくお願いします。

  • 窓関数の処理の仕方によるFFT結果の違いについて

    ある1周期分のデータX(n)={X0,X1,・・・,Xn-1}があるとき(nは2の乗数)、ハミング窓関数K=0.54-0.46cos(2πn/N-1)を以下のように処理 (1)Y(n)=K*X(n) (2)Y(n)=K*X(n)+A (3)Y(n)=K*{X(n)+A} ただし、AはX(n)の平均値。 ここで、Y(n)をFFTした場合の違いについてシミュレーションしましたが、振幅スペクトルの結果が次のようになりました。  (1)と(2)は0次のみ異なる(1次以降は同じ)  (2)と(3)((1)と(3))は0次と1次が異なる(2次以降は同じ) 【質問1】 (2)と(3)の2次以降が同じになる理由はどうしてでしょうか。

  • 離散フーリエ変換のプログラムについて

    今DFTのプログラムをC言語で書いているんですが、うまく動いてくれません。 DFTの式はX(k)=1/N{Σx(k)*e^(-j2πkn/N)}のシグマの中をn=0からN-1まで足し合わせればいいと思っているのですがちがいますか。 下にプログラムを書きますのでお願いします。 void DFTcore(double x[],int N,double A[],double fai[],double a_rl[],double a_im[]){     x[]はデータ     Nはデータ数     A[]は振幅     fai[]は位相差     a_rl[]は実数部、a_im[]は虚数部 double w,a; int k,n,t; for(k=0;k<N;k++){  //初期化    a_rl[k]=0.0;  //実数部    a_im[k]=0.0;  //虚数部   }     時間間隔は1秒 for(k=0;k<N;k++){       for(n=0;n<N;n++){   w=((2*PI)/N)*k*n;         a_rl[k]=a_rl[k]+x[k]*cos(w);   a_im[k]=-a_im[k]+x[k]*sin(w);  }   a_rl[k]=a_rl[k]/(double)N;実数部   a_im[k]=-a_im[k]/(double)N; 虚数   A[k]=sqrt(a_rl[k]*a_rl[k]+a_im[k]*a_im[k]); //振幅    fai[k]=atan(a_im[k]/a_rl[k]); //位相 } }

  • データの分析 118

    nを2以上の自然数とする。次の問いに答えよ。 (1)変量xの値がx(1),x(2),・・・,x(n)であるとし f(a)=1/nΣ[k=1,n](x(k)-a)^2 とする。f(a)を最小にするaはx(1),x(2),・・・,x(n)の平均値で、そのときの最小値はx(1),x(2),・・・,x(n)の分散であることを示せ。 (2)cを定数として、変量y,zのk番目のデータの値が y(k)=k(k=1,2,・・・,n),z(k)=ck(k=1,2,・・・,n) であるとする。このときy(1),y(2),・・・,y(n)の分散がz(1),z(2),・・・,z(n)の分散より大きくなるためのcの必要十分条件を求めよ。 (3)変量xのデータの値がx(1),x(2),・・・,x(n)であるとし、その平均値をxバーとする。新たにデータを得たとし、その値をx(n+1)とする。x(1),x(2),・・・,x(n),x(n+1)の平均値をx(n+1),xバーおよびnを用いて表せ。 (4)次の40個のデータと平均値、分散、中央値を計算すると、それぞれ、ちょうど40,670,35であった。 120,10,60,70,30,20,20,30,20,60 40,50,40,10,30,40,40,30,20,70 100,20,20,40,40,60,70,20,50,10 30,10,50,80,10,30,70,10,60,10 新たにデータを得たとし、その値が40であった。この時、41個のすべてのデータの平均値、分散、中央値を求めよ。ただし、得られた値が整数でない場合は、少数第1位を四捨五入せよ。 この問題を解いてください。お願いします。

  • DCモータの騒音分析について

    ブラシ付きDCモータの騒音分析をしていますが、少々わからないことがあって質問いたします。 普通騒音計を使用し、A特性を用いてFFTアナライザで騒音測定しました。 出力をx-y軸それぞれ周波数-騒音レベル[dB]としたところ、 ある周波数k[Hz]で大きなピークがありました。 このピークの原因は、モータのブラシとコミューテータの接触によるものだと考えるのが普通だと思い、この接触で発生する騒音周波数を算出しました。 この騒音周波数nとさきほどの騒音ピーク時の周波数kの関係は、 k = 2 * n となりました。 よく、基音周波数の整数倍にピークがくることがあるということを聞きますが、これはどのような理論で導き出されているのでしょうか? やはり共振・共鳴とかの話でしょうか?しかし、今回の場合、ブラシとコミューテータの接触する時の周波数から導き出した値の倍数なので、 私の中では、共振・共鳴とは考えずらいのですが。 ブラシとコミューテータの接触する音(時系列では矩形波?)をFFTで周波数領域に落とし込むと、整数倍の周波数成分帯に反応する原理があるのでしょうか?もしこれだとしたら、詳しく教えていただきたいです。よろしくお願いします。

  • 離散フーリエ変換について

    10 REM dft 20 N=8 30 DIM A(N),B(N),X(N) 40 FOR I=0 TO N-1 50 READ X(1) 60 NEXT I 70 P=6.283/N 80 FOR K=0 TO N-1 90 A(K)=0:B(K)=0 100 FOR J=0 TO N-1 110 A(K) =A(K)+X(j)*COS(P*J*K) 120 B(K) =B(K)-X(J)*SIN(P*J*K) 130 NEXT J 140 NEXT K 150 FOR I=0 TO N-1 160 Y=SQL(A(I)^2+B(I)^2) 170 LPRINT I;:LPRINT USING "###,###";A(I),B(I),Y 180 NEXT I 190 DATA 1,1,1,1,0,0,0,0 このDFTプログラムをC言語に直したいのですがよく分かりません; お願いします@@;

  • 面積についての矛盾(?)

    はじめまして。 現在高2です。 他の掲示板で質問をしていたのですが、満足のいく返信が来なかったためここで質問させてください。 0≦x≦1で、f(x)≧0を満たす関数f(x)について、f(x),x軸,y軸,x=1で囲まれた範囲の面積Sは、 S:=∫_0^1 f(x)dx=lim_{n→∞}Σ_{k=1}^n {1/n f(k/n)} と定義できますが(多分あってると思います)、 a≠b,b→a として、y=b,x軸,y軸,x=1で囲まれた範囲の面積Sを考えれば、 明らかにS=b ところで、上の定義に従えば、 S=lim_{n→∞}Σ_{k=1}^n {b/n} で、これは n→∞について、 Σ_{k=1}^n {b/n}→bが成り立つことを示している。 ここで、b→aなので、 Σ_{k=1}^n {b/n}→b→a Σ_{k=1}^n {b/n}→a つまり、 n→∞について、 Σ_{k=1}^n {b/n}→a が成り立つことを示している。 よって S=lim_{n→∞}Σ_{k=1}^n {b/n}=a ゆえに、S=b=a これは仮定a≠bに反する。 いったいどこがおかしいのでしょうか。。 数式が見にくくてすみません。

  • 【エクセル】株価をフーリエ変換し、波を三つ作りたい

    ※実際はドル円ですが、タイトルの文字数制限の関係で株価と表記することにします エクセルで株価の終値のデータ、128個用意したものをFFTしました。 A0=FFTの実部÷128 // 解析する周期2L=128 A1~A127=FFTの実部÷(128÷2) // 解析する周期2L÷2=128÷2 B1=(-1)*FFTの虚部÷(128÷2) // 解析する周期2L÷2=128÷2 振幅A゜1~A゜63(ナイキスト周波数を考慮して前半A゜0~A゜63のうち、A゜0を除いた)についての棒グラフを作成(→振幅スペクトルのグラフ) ※振幅A゜nはフーリエ係数Anと区別するために別の記号を用いた 振幅スペクトルで大きいほうから3つ(取り出し方はVLOOKUPで取り出したとする)を取り出し、それぞれの波を株価グラフ(折れ線)上に描きたいのですが、どうすれば良いのでしょうか? そこだけが分かりません。 画像が荒くてすみません。