2

matlabヘルプページにあるのとまったく同じように対数チャープを作成します。

t = 0:0.001:10;      % 10 seconds @ 1kHz sample rate
fo = 10; f1 = 400;   % Start at 10Hz, go up to 400Hz
X = chirp(t,fo,10,f1,'logarithmic');
figure(2);
spectrogram(X,256,200,256,1000,'yaxis');

スペクトログラム

次に、他のアプリケーションで機能する次のコードを使用して、周波数領域に移動します。

fft_prep = fftshift(fft(X));
fft_mag = abs(fft_prep);
pos_fft = fft_mag(1:ceil(length(fft_mag)/2));
db_fft = 20*log10(pos_fft);
figure(1);
plot(db_fft);

そして、次のグラフが1kHz〜5kHzで刺激的であるように見えるのを見て驚いた。

SpectrumAnalyzer

私はmatlabのチャープ機能に精通しておらず、誰かが私が見逃している明らかな何かを見たのではないかと思っていました。他のポインタは大歓迎です。

4

4 に答える 4

7

機能に問題はありませんchirp...

ベクトルインデックスではなく、周波数値に対して db_fft をプロットするだけです =)。

plot(linspace(fo,f1,length(db_fft)), db_fft);

ここに画像の説明を入力

また、他の FFT メソッドを使用して信号の FFT を計算することもテストしましたが、それらも 0 ～ 400 Hz の範囲を示しています。

アップデート：

IMO、dBまたはパワー（ピリオドグラム）でプロットしない方が視覚的に簡単だと思います。以下は、時間領域信号の FFT を計算するための優れた例と私の goto メソッドです: mathworks.se/help/matlab/ref/fft.html

応答：

少し考えた後、私は上記の私の答えが間違っていることに同意しますが、あなたが言う理由ではありません. 周波数領域の x 軸は、チャープの実際の長さ (または半分、2 倍、またはそのようなもの) に移動しないでください。周波数領域の x 軸は、信号のサンプリング レートの半分 (Fs/2) に移動する必要があり、希望/希望する最大周波数の 2 倍のサンプリング周波数で信号をサンプリングしていることを確認する義務があります。解決。

言い換えれば、FFT が時間領域信号の長さと同じ/2 倍/半分であると仮定するのは正しくありません。これは、FFT を表すために任意の数の周波数ビンを選択できるためです。ベスト プラクティスは長さ = N^2 ( 2 の累乗) を使用して計算を高速化します。FFT を計算するときに、なぜ時間値を知る必要があるのでしょうか。あなたはそうしない！サンプリング周波数のみが必要です (Fs = 0.001 ではなく、Fs = 1000 btw に設定する必要があります)。

上記の私の答えは正しくありません。次のようにする必要があります。

plot(linspace(0, Fs/2, length(db_fft)), db_fft)

Fs/2 の代わりに、length(t)/(2*Tfinal) と書きました。これは (ほぼ) Fs/2 と同じ値ですが、正しい方法ではありません =)。

これが私のgoto FFTメソッドです（値はdBではありません）。

function [X,f] = myfft(x,Fs,norm)
    % usage: [X, f] = myfft(x,Fs,norm);
    %        figure(); plot(f,X)
    % norm: 'true' normalizes max(amplitude(fft))=1, default=false.
    if nargin==2
        norm=false;
    end
    L = length(x); NFFT = 2^nextpow2(L);
    f = Fs/2*linspace(0,1,NFFT/2+1);
    %f =0:(Fs/NFFT):Fs/2;
    X = fft(x,NFFT)/L; X = 2*abs(X(1:NFFT/2+1));
    if norm==true; X = X/max(abs(X)); end
end

[Xfft, f] = myfft(X,Fs); の結果のグラフは次のとおりです。プロット (f、Xfft); NyQuist の定理に従って、リターン周波数ビン ベクトルは max(f) = Fs/2 であることに注意してください (Fs/2 よりも高い周波数は解決できません)。

ここに画像の説明を入力

于 2013-03-14T07:22:54.217 に答える
0

いくつかのエラーがありましたが、すべてが修正されたわけではありません。コードをもう少し試した後、ここに私が思いついた解決策があります。

セットアップをよりモジュール化するために、さらにいくつかの変数を追加しました。

Tfinal = 10;
Fs = 1/1000;
t = 0:Fs:Tfinal;      % 10 seconds @ 1kHz sample rate
fo = 10; f1 = 400;   % Start at 10Hz, go up to 400Hz
X = chirp(t,fo,Tfinal,f1,'linear');

マグニチュードとリンスペースをプロットすると、リンスペースは、低周波から高周波までだけでなく、実際のチャープ信号の長さに一致する必要があります。ベクトルtの長さは 1000 であり、FFT の後で正の半分を使用するため、チャープ信号の長さは 500 であり、リンスペースf1が示唆する 400 ではありません。

fft_prep = fftshift(fft(X));
fft_mag = abs(fft_prep);
pos_fft = fft_mag(ceil(length(fft_mag)/2)+1:length(fft_mag));
db_fft=20*log10(pos_fft);
figure(1);
plot(linspace(0,length(t)/(2*Tfinal),length(db_fft)), db_fft);

また、前述の修正を行って、FFT のマイナス側ではなくプラス側を取得しました。これはプロットします：

ここに画像の説明を入力

これは、10Hz-400Hz を励振するチャープを正確にグラフ化したものです。極端な場合に行って線形にすることで、より明確に見ることができます。linspace 補正なしで、サンプリング周波数を 100、範囲を 10 ～ 25 にしてみます。

ここに画像の説明を入力

修正後：

ここに画像の説明を入力

于 2013-03-15T19:02:50.683 に答える
0

ちなみに、元のチャープ振幅を FFT から抽出すると便利な場合があります。たとえば、あるデバイスの周波数でオーディオ スイープを行っていて、各周波数での応答の振幅と位相を知りたい場合 (Pspice のように)。つまり、振幅は a^2 = abs(FFT)^2 *4 * チャープの帯域幅 /(Fs * N) で、Fs はサンプル周波数、N は FFT のポイント数です。たとえば、200 ～ 400Hz のチャープの帯域幅は 200Hz です。

導出を知りたい場合は、Parseval の定理から始めます: 時間信号の平均 2 乗 = PSD の下の面積。したがって、a^2/2 = sum(abs(FFT)^2) / N^2 ここで、a は掃引周波数信号 (チャープなど) のピーク振幅です。チャープが周波数分布において対数ではなく線形である場合、FFT は上記のプロットのようにフラットになります。したがって、合計を Nb * abs(FFT)^2 / N^2 で置き換えることができます。ここで、Nb はチャープが占有する周波数ビンの数であり、abs(FFT) は FFT の大きさであり、すべての値で同じです。チャープが占める周波数ビン。1 つの周波数ビンの帯域幅が Fs/N であるという事実を使用すると、チャープの帯域幅は Nb * Fs /N であることがわかります。上記の結果は、ここから簡単に導出できます。

于 2014-05-16T02:58:49.340 に答える