728x90

 0. [matlab] - Continuous-Time Fourier Transformation





 1. 풀이


이전에 코드들 (https://kyunstudio.tistory.com/281(퓨리에 시리즈), https://kyunstudio.tistory.com/278(적분))를 활용해서 구현을 해준 문제였다.

실제 코드를 실행해보면 재밌는 그래프들을 얻을 수 있을 것이라 생각한다.


기본 함수를 즉, origin function을 function handle로 만들어 주는 것이 가장! 가장! 중요한 포인트이다.

제대로 구현만 완성하면 어떤 함수에 대해서도 fourier transform 테스트가 가능하다.



 2. 소스코드


1.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
clear;
% (i)
T_1 = 10;
f_t = @(t) rectpuls(t, T_1);        %rectpuls함수 활용해 사각파를 만든다.
 
% 방법 1. C_k를 활용해 F(w)를 구한다.
% 이때 T_0 = T_1*1000을 취해준다.(극한을 취해주어야하지만, 1000배를 해주었다.)
T_0 = T_1*1000;
dw = 2*pi/T_0;
= -2000:1:2000;     % k= -inf:1:inf 이어야 하지만, 계산상의 문제로 4000개만 계산하였다.
C_1 = my_fourier_coefficient(f_t, T_0, k);
F_w_1 = C_1*2*pi/dw;
= 2*pi*k/T_0;         % w는 k에 의해 결정
 
 
%%%%% 그래프 출력
figure(1); subplot(211);stem(w, abs(F_w_1));
title('Amplitude')
xlabel('w');
ylabel('|F(w)|');
subplot(212); stem(w, angle(C_1)*180/pi);
title('Angle')
xlabel('w');
ylabel('degree of F(w)');
 
 
 
% 방법 2. F(w)를 적분을 활용해 직접 만든다.
f_w = @(w) my_fourier_transformation(f_t, w, T_1);  % w에 대한 함수를 만들어준다.
= -20:0.01:20;    % 이때, w는 연속적으로 변화하는 값으로 설정
                    % 이러한 값은 T_0가 inf인 경우이고
                    % 이산적인 결과가 아닌 연속적인 그래프가 나오게 된다.
                    
                    
= -T_1*10:T_1/100:T_1*10;
= f_t(t);
y_2 = f_w(w);
 
%%%%% 그래프 출력
figure(2);subplot(311); plot(t, y);
title('origin')
xlabel('t');
ylabel('f(t)');
 
subplot(312); plot(w, abs(y_2));    % 그래프 출력에 stem이 아닌 plot를 활용
title('fourier transformation(amplitude)')
xlabel('w');
ylabel('|F(w)|');
 
subplot(313); plot(w, angle(y_2)*180/pi);
title('fourier transformation(Angle)');
xlabel('w');
ylabel('degree of F(w)f(w)');
 
 
function a = myintegral(fun, gap, min, max)
    a = 0;      % 결과로 반환될 a를 0으로 초기화
    for i = min : gap : max
        a = a + fun(i)*gap;      % 반복문 활용한 수치 적분
    end
end
 
function ret = my_fourier_transformation(fun, w, T)   
    FUN = @(t) fun(t).*exp(-i*w.*t);     % 2*pi*k/T_0를 w로 치환
    ret = myintegral(FUN, T / 10000,-T,T);    % 직접 만든 적분을 활용해 F(w)을 구한다.
end
 
function ret = my_fourier_coefficient(fun, T_0, k)   
    FUN = @(t) fun(t).*exp(-i*k*(2*pi/T_0).*t);             % 퓨리에 변환을 활용한다.
    ret = (1/T_0)*myintegral(FUN, 0.1-T_0/2, T_0/2);     % 직접 만든 적분을 활용해 C_k을 구한다.
end
cs



2.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
clear;
% (i)
T_1 = 100;
= -T_1*10:T_1/100:T_1*10;
f_t = @(t) rectpuls(t, T_1);        %rectpuls함수 활용해 사각파를 만든다.
 
% 방법 1. C_k를 활용해 F(w)를 구한다.
% 이때 T_0 = T_1*1000을 취해준다.(극한을 취해주어야하지만, 1000배를 해주었다.)
T_0 = T_1*1000;
dw = 2*pi/T_0;
= -2000:1:2000;     % k= -inf:1:inf 이어야 하지만, 계산상의 문제로 4000개만 계산하였다.
C_1 = my_fourier_coefficient(f_t, T_0, k);
F_w_1 = C_1*2*pi/dw;
= 2*pi*k/T_0;         % w는 k에 의해 결정
 
 
%%%%% 그래프 출력
figure(1); subplot(211);stem(w, abs(F_w_1));
title('Amplitude')
xlabel('w');
ylabel('|F(w)|');
subplot(212); stem(w, angle(F_w_1)*180/pi);
title('Angle')
xlabel('w');
ylabel('degree of F(w)');
 
 
 
% 방법 2. F(w)를 적분을 활용해 직접 만든다.
F_w = @(w) my_fourier_transformation(f_t, w, T_1);  % w에 대한 함수를 만들어준다.
= -20:0.01:20;    % 이때, w는 연속적으로 변화하는 값으로 설정
                    % 이러한 값은 T_0가 inf인 경우이고
                    % 이산적인 결과가 아닌 연속적인 그래프가 나오게 된다.
= f_t(t);
y_2 = F_w(w);
 
%%%%% 그래프 출력
figure(2);subplot(311); plot(t, y);
title('origin')
xlabel('t');
ylabel('f(t)');
 
subplot(312); plot(w, abs(y_2));    % 그래프 출력에 stem이 아닌 plot를 활용
title('fourier transformation(amplitude)')
xlabel('w');
ylabel('|F(w)|');
 
subplot(313); plot(w, angle(y_2)*180/pi);
title('fourier transformation(Angle)')
xlabel('w');
ylabel('degree of F(w)f(w)');
 
 
 
function a = myintegral(fun, gap, min, max)
    a = 0;      % 결과로 반환될 a를 0으로 초기화
    for i = min : gap : max
        a = a + fun(i)*gap;      % 반복문 활용한 수치 적분
    end
end
 
function ret = my_fourier_transformation(fun, w, T)   
    FUN = @(t) fun(t).*exp(-i*w.*t);     % 2*pi*k/T_0를 w로 치환
    ret = myintegral(FUN, T / 10000,-T,T);    % 직접 만든 적분을 활용해 F(w)을 구한다.
end
 
function ret = my_fourier_coefficient(fun, T_0, k)   
    FUN = @(t) fun(t).*exp(-i*k*(2*pi/T_0).*t);             % 퓨리에 변환을 활용한다.
    ret = (1/T_0)*myintegral(FUN, 1-T_0/2, T_0/2);     % 직접 만든 적분을 활용해 C_k을 구한다.
end
cs





 3. 참고




질문이나 지적 있으시면 댓글로 남겨주세요~

도움 되셨으면 하트 꾹!


+ Recent posts