푸리에 급수의 최대급수 차수를 (N=2,5,10,20)으로 변경하면서 x(t) 신호와 푸리에 급수로 표현되는 신호를 Matlab 을 이용하여 하나의 Graph 에서 비교하시오.
- Matlab 으로 Graph 도시할 때, [-T, 2T] 구간에서 그래프를 도시할 것
- T1 = 0.5msec, T = 1msec, A = 1
close all; clear all;
A = 1;
T = 1e-3;
w0 = 2 * pi / T;
degree = [2 5 10 20]; % 최대차수 설정
d1 = 0;
d2 = 0;
d3 = 0;
d4 = 0;
a0 = A / 2;
t = linspace(-1e-3 , 2e-3, 10000);
for n = 1 : degree(1)
% a, b는 푸리에 계수
a(n) = 2 * A * sin(n*w0*T/2) / (T * n * w0);
b(n) = 2 * A * (1 - cos(n*w0*T/2)) / (T * n * w0);
ct = a(n) * cos(n*w0*t) + b(n) * sin(n*w0*t);
% ct는 푸리에 계수를 포함한 푸리에 급수 식(아직 a0X)
d1 = d1 + ct;
end
for n = 1 : degree(2)
a(n) = 2 * A * sin(n*w0*T/2) / (T * n * w0);
b(n) = 2 * A * (1 - cos(n*w0*T/2)) / (T * n * w0);
ct = a(n) * cos(n*w0*t) + b(n) * sin(n*w0*t);
d2 = d2 + ct;
end
for n = 1 : degree(3)
a(n) = 2 * A * sin(n*w0*T/2) / (T * n * w0);
b(n) = 2 * A * (1 - cos(n*w0*T/2)) / (T * n * w0);
ct = a(n) * cos(n*w0*t) + b(n) * sin(n*w0*t);
d3 = d3 + ct;
end
for n = 1 : degree(4)
a(n) = 2 * A * sin(n*w0*T/2) / (T * n * w0);
b(n) = 2 * A * (1 - cos(n*w0*T/2)) / (T * n * w0);
ct = a(n) * cos(n*w0*t) + b(n) * sin(n*w0*t);
d4 = d4 + ct;
end
ft1 = a0 + d1;
ft2 = a0 + d2;
ft3 = a0 + d3;
ft4 = a0 + d4;
%d1, d2, d3, d4에 a0값을 더해 푸리에 급수식을 완성
plot(t, ft1, t, ft2, t, ft3, t, ft4);
legend('N = 2', 'N = 5', 'N = 10', 'N = 20');
xlabel('time(t)');
ylabel('Fourier');
푸리에 급수로 표현을 할 때 최대 차수 설정을 바꿔가면서 파형을 도출하는 실습이었다. 사각펄스 함수를 푸리에 급수로 표현을 하고 차수를 늘려감에 따라 비교분석을 해보았다. 차수가 낮을 때는 표현되는 주파수 성분이 적어 원래의 함수 파형과는 많이 다른 모습이었다.
하지만 최대차수를 늘리면 푸리에 급수에 의해 더 많은 삼각함수 성분(기저함수)들이 증가하고 이에 포함되는 주파수 성분도 증가하여 원래 함수와 점점 가까워지는 것을 확인할 수 있었다.
Question.
위의 b)에서 푸리에 급수의 최대급수 차수를 (N=2,5,10,20)으로 변경하는 경우에 원신호 x(t)와의 RMS 오차 를 비교하는 Matlab 프로그램을 구성하고, 그 값을 제시하시오.
close all; clear all;
A = 1;
T = 1e-3;
T1 = T /2;
w0 = 2 * pi / T;
degree = 20;
d = 0;
t = linspace(-1e-3 , 2e-3, 9000); %샘플 개수를 9000
a0 = A / 2;
for n = 1 : degree
a(n) = 2 * A * sin(n*w0*T/2) / (T * n * w0);
b(n) = 2 * A * (1 - cos(n*w0*T/2)) / (T * n * w0);
ct = a(n) * cos(n*w0*t) + b(n) * sin(n*w0*t);
d = d + ct;
end
ft = a0 + d;
%한 주기동안 3000개의 샘플값을 가지므로, 그의 절반은 1500
t1 = linspace(-T, -T1, 1500);
t2 = linspace(-T1, 0, 1500);
t3 = linspace(0, T1, 1500);
t4 = linspace(T1, T, 1500);
t5 = linspace(T, 3 * T / 2, 1500);
t6 = linspace( 3 * T / 2 , 2 * T , 1500);
tk = [t1 t2 t3 t4 t5 t6]; %단위계단함수의 샘플된 t
x1 = ones(size(t1));
x2 = zeros(size(t2));
x3 = ones(size(t3));
x4 = zeros(size(t4));
x5 = ones(size(t5));
x6 = zeros(size(t6));
%ones, zeros 함수로 단위 계단 함수 표현
x = [x1 x2 x3 x4 x5 x6];
for i = 1 : 3000
d(i) = (ft(i) - x(i))^2;
end
add = sum(d);
rmse = sqrt(add/3000);
trmse = 3 * rmse
% RMSE를 구하기 위한 코드
RMSE(평균 제곱근 편차) : 추정 값 또는 모델이 예측한 값과 실제 환경에서 관찰되는 값의 차이를 다룰 때 흔히 사용하는 측도이다. 사각펄스와 푸리에 급수로 표현된 펄스의 평균 제곱근 편차를 구하는 실습이었다. 최대 차수를 늘릴수록 기존 파형과 비슷한 파형이 나오기 때문에 차수를 늘릴수록 오차도 줄어야할 것이다. 사각펄스를 만들기 위해 ones와 zeros를 사용하여 배열을 만들고 합쳤다. 3주기만큼의 펄스를 만들었기 때문에 1주기의 RMSE를 만들고 3을 곱해주었다. 사각펄스(x)와 푸리에급수로 표현된 펄스(ft)를 선언하고 두 값 차이의 제곱을 d(i)로 저장한 후에 그 값을 sum하고 한 주기의 sample값인 3000으로 나눠주었다. 최종적으로 이 값에 루트를 씌워 RMSE값을 도출하였다. 실제로 최대차수를 2, 5, 10, 20 으로 늘려감에 따라 RMSE 값이 감소하는 것을 확인할 수 있었다.
