0

フーリエ変換するために、基本周期2秒の矩形パルスxを作成しました

t=-2:0.01:2; 
dt=t(2)-t(1); %increment of time
fs=1/dt % sampling rate 
n=length(t) %number of samples
X=fftshift(fft(x,n))/n %fourier transform of x%

f=linspace(-fs/2,fs/2, n) %making frequency axis

X_angle=angle(X); %the phase of X

私は、位相スペクトルが交互に -pi/2 と pi/2 になることを期待していました。

しかし、グラフ（評判が悪いために投稿できないのは残念です）

X_angle は、周波数が増加するにつれて徐々に増加し、-pi から pi の範囲になることを示しています

マグニチュードスペクトルをプロットするとうまくいきました

plot(f,X_mag), X_mag=abs(X)

私はいくつかの大きさの調整をしなければならないと思うangle(X)

しかし、角度は X の大きさと無関係ではありませんか?

絶対周波数が高くなると X_angle の絶対値が大きくなる理由がわかりません。

4

1 に答える 1

1

コードでいくつかの間違いを犯しました。

  1. コードで x の値を出力しませんでした。対称矩形パルスを分析したいと思います。ただし、サンプル数は奇数です。これは、信号が対称ではないことを意味し、理論的な教科書の値と比較して小さな違いが生じます。

    1. 方形パルスの帯域幅は無限です。ただし、FFT の場合、サンプリング レートは信号の最高周波数の少なくとも 2 倍である必要があります。これは、サンプリング レートが少なくとも 2*無限大でなければならないことを意味します。残念ながら、それはできません。誰もそれを行うことはできません。その結果、エイリアシングが発生します。つまり、結果にエラーが含まれます。幸いなことに、矩形関数の場合、この影響はsinc関数を使用して補正できます。

    2. 何かを正しく行えば、正しい教科書係数が得られます。入力信号が (N=10) 11111-1-1-1-1-1 のような形式の場合、関数は奇関数です。これは、f(-t) = -f(t) を意味します。この場合、四角形は一連の正弦関数によって作成できます。理論上の関数は次のとおりです。

    f(t) = 4/pi( sin(wt) + 1/3 sin(3wt) + 1/5 sin(5wt) + 1/7 sin(7wt) ...)

    sin(2wt) や sin(4wt) のような偶数の周波数はありません。これは、スペクトル内の 2 つおきの周波数の値がゼロであることを意味します。数値ノイズが原因で、これらの値は正確にゼロではなく、ゼロに近くなっています。これらの値から位相を計算すると、無意味な値が生成されます。他の周波数は、正弦関数のフーリエ変換値です。正弦関数の FFT は純粋な虚数です。合計のすべての要素が同じ符号を持つため、すべての角度は同じ値 (pi/2) になります。

以下に、変更されたコードを示します。

close all;
clear all;
clc;

t=-2:0.01:(2-0.01); 
dt=t(2)-t(1); %increment of time
fs=1/dt; % sampling rate 
n=length(t); %number of samples
x = [ones(1,n/2), -ones(1,n/2)];
X=fft(x)/n; %fourier transform of x%

f=0:n-1; %making frequency axis
sincComp = @(N) (exp(-i*pi*(f/N)).*sinc(f/N)).';
%  Perform the compensation
comp = transpose(sincComp(n));
Xcorrected = X.*comp;

X_angle=angle(Xcorrected(2:2:n)); %the phase of X

figure;
subplot(3,1,1);
hold on;
stem(f(1:10), pi*abs(X(1:10))/2, 'ob', 'LineWidth', 3);
plot(f(1:10), pi*abs(Xcorrected(1:10))/2, 'or', 'LineWidth', 3);
hold off;
grid on;
title('Absolute value of FFT result', 'FontSize', 18);
xlabel('frequency', 'FontSize', 18);
ylabel('abs', 'FontSize', 18);
legend(['FFT'], ['FFT compensated']);
subplot(3,1,2);
hold on;
stem(f(1:10), pi*real(X(1:10))/2, 'ob', 'LineWidth', 3);
plot(f(1:10), pi*real(Xcorrected(1:10))/2, 'or', 'LineWidth', 3);
hold off;
grid on;
title('Real value of FFT result', 'FontSize', 18);
xlabel('frequency', 'FontSize', 18);
ylabel('real', 'FontSize', 18);
subplot(3,1,3);
hold on;
stem(f(1:10), pi*imag(X(1:10))/2, 'ob', 'LineWidth', 3);
plot(f(1:10), pi*imag(Xcorrected(1:10))/2, 'or', 'LineWidth', 3);
hold off;
grid on;
title('Imaginary value of FFT result', 'FontSize', 18);
xlabel('frequency', 'FontSize', 18);
ylabel('imag', 'FontSize', 18);

FFT と亜鉛補償 FFT の結果

FFT と亜鉛補償 FFT の結果

最初の 10 個の非ゼロ周波数の位相値:

>> X_angle(1:10)*2/pi
ans =
-1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000
于 2014-09-30T20:37:32.643 に答える