어떻게 정확하게 matlab에 푸리에 시리즈의 위상 스펙트럼을 그릴 수 있습니까?

Question

가 나는 위상 스펙트럼 대체물을 기대 하였다

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 

푸리에 변환에, 기본주기 2초과 직사각형 펄스 X를 만들었다 -pi/2 PI/2

하지만 그래프 (내가 인해 내 평판의 부족)

이 X_angle 점차 주파수가 증가함에 따라 증가 나에게 보여줍니다에 게시 할 수 없습니다 나쁘지 내가 magnitu를 꾸몄다 때

는 그냥 괜찮 았는데 -pi에서 파이의 범위 드 스펙트럼은,

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

와 나는 angle(X)

일부 크기 조정을해야하지만, 각도가 X의 크기와 무관하지 않다 생각하십니까?

왜 나는 절대 주파수가 증가함에 따라 X_angle의 절대 값이 증가하는지 알지 못합니다.

Answer 1

코드에서 몇 가지 실수를 저질렀습니다.

1. 코드에서 x의 값을 인쇄하지 않았습니다. 대칭 직사각형 펄스를 분석한다고 가정합니다. 그러나 샘플 수는 홀수 값입니다. 이것은 신호가 대칭이 아니므로 의 차이가 이론상의 텍스트 북 값과 비교하여 작음을 의미합니다.

   2. 직사각형 펄스의 대역폭은 무한합니다. 그러나 FFT의 경우 샘플링 속도는 신호의 최고 주파수 인 의 두 배 이상이어야합니다. 즉, 샘플링 속도는 최소 2 * 무한대 여야합니다. 불행, 할 수 없다. 아무도 그렇게 할 수 없습니다. 결과로 앨리어싱이 발생합니다. 즉, 결과에 오류가 있음을 의미합니다. 좋은 소식은 직사각형 함수의 경우이 효과가 sinc 함수의 도움으로 보상 될 수 있다는 것입니다.

   3. 올바른 것을 만들면 올바른 교과서 계수를 얻을 수 있습니다. 입력 신호의 형식이 (N = 10) 11111-1-1-1-1-1 인 경우 함수는 홀수 함수입니다. 이것은 f (-t) = -f (t)를 의미합니다. 이 경우 직사각형은 일련의 사인 함수로 구성 할 수 있습니다. 이론적 함수이다 :

   F (t) = 4/PI (SIN (중량) + 1/3 죄 (3 중량) + 1/5 죄 (5 중량) + 1/7 죄 (7 중량). ..)

   sin (2wt) 또는 sin (4wt)와 같은 짝수의 빈도가 없습니다. 이것은 스펙트럼의 모든 두 번째 주파수가 0이라는 것을 의미합니다. 수치 노이즈로 인해이 값은 정확히 0은 아니지만 0에 가깝습니다. 이 값에서 위상을 계산하면 의미없는 값이 생성됩니다. 다른 주파수는 사인 함수의 푸리에 변환 값입니다. 사인 함수의 FFT는 순수한 허수입니다. 합계의 모든 요소에는 동일한 부호가 이므로 모든 각도는 같은 값을 갖습니다. 즉 pi/2입니다. 당신이 수정 된 코드를 찾아 아래

: FFT의

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를

처음 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