[matlab] Continuous-time Fourier Series (연속 시간, 퓨리에 시리즈, function handle 활용, 함수의 신기한 활용,신호 및 시스템, Signals and Systems)

2020. 7. 19. 23:25

728x90

0. [matlab] - Continuous-time Fourier Series

1. 풀이

이번 문제에서는 익명 함수를 함수에 바로 대입이 가능하도록 만들어 준 것이 핵심 포인트이다.

문제를 해결하며, 어떠한 함수에 대해서도 퓨리에 시리즈를 얻을 수 있도록 익명 함수(function handle)을 선언하여 이 익명 함수를 함수의 인자로 전달하는 것이 핵심 포인트였다.

또한 함수의 return 또한 구간 크기만큼 반환을 하도록 만들어주었다.

원하는 구현을 하도록 오랜 시간이 걸렸지만, 만족했던 코딩 :)

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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
clear;
% (i) Amplitude plot를 출력
 
% function handle 선언
x_1 = @(t) square_wave(t, 1, 20);
x_2 = @(t) sawtooth_wave(t, 1, 0.02);
 
k = -10:1:10;
 
% 1. input : square wave
C_1 = my_fourier_coefficient(x_1, 20, k);       % Ck를 구한다.(k = -10:1:10)
t = 0:100/10000:100;
x_1_1 = zeros(1,10001);                         % 10001개의 x_1_1배열을 생성
for a = 1:1:10001
    x_1_1(a) = x_1(t(a));                       % x_1_1에 squere wave를 담아준다.
end
 
 
% 그래프 출력
figure(1); subplot(211);plot(t, x_1_1);
title('Square wave')
xlabel('t');
ylabel('x(t)');
figure(1); subplot(212);stem(k, abs(C_1));
title('Amplitude')
xlabel('k');
ylabel('|C_k[n]|');
 
 
 
% 2. input : sawtooth wave
C_2 = my_fourier_coefficient(x_2, 0.02, k);     % Ck를 구한다.(k = -10:1:10)
t = 0:0.1/10000:0.1;
x_2_2 = x_2(t);                                 % x_2_2에 sawtooth wave를 담아준다.
 
% 그래프 출력
figure(2); subplot(211);plot(t, x_2_2);
title('Sawtooth wave')
xlabel('t');
ylabel('x(t)');
figure(2); subplot(212); stem(k, abs(C_2));
title('Amplitude')
xlabel('k');
ylabel('|C_k[n]|');
 
 
 
function a = myintegral(fun, gap, min, max)
    % function handle 검사
    if ~isa(fun,'function_handle')
        error(message('MATLAB:integral:funArgNotHandle'));  % 오류메세지 출력 후 프로세스 종료
    end
    
    a = 0;                      % 결과로 반환될 a를 0으로 초기화
    
    for i = min : gap : max
        a = a + fun(i)*gap;     % 반복문 활용한 수치 적분
    end
end
 
% 결과를 배열에 담아 return 하는 convolution 함수
% x,h를 배열로 입력
function ret = my_fourier_coefficient(fun, T_0, k)
    % function handle 검사
    if ~isa(fun,'function_handle')
        error(message('MATLAB:integral:funArgNotHandle'));  % 오류메세지 출력 후 프로세스 종료
    end
    
 
    FUN = @(t) fun(t).*exp(-i*k*(2*pi/T_0).*t);             % 퓨리에 변환을 활용한다.
    ret = (1/T_0)*myintegral(FUN, T_0 / 10000, 0, T_0);     % 직접 만든 적분을 활용해 C_k을 구한다.
        
end
 
function ret = my_fourier_series(fun, N, T_0, t)
    % function handle 검사
    if ~isa(fun,'function_handle')
        error(message('MATLAB:integral:funArgNotHandle'));  % 오류메세지 출력 후 프로세스 종료
    end
    
    ret = 0;
    
    for k=-N:1:N
        C_k = my_fourier_coefficient(fun, T_0, k);
        ret = ret + C_k*exp(i*k*(2*pi/T_0)*t);
    end
end
 
function ret = square_wave(t,X_0, T_0)
    T = mod(t,T_0);     % T는 주기로 나눈 나머지
    if(T < (T_0 / 2))
        ret = X_0;
    else
        ret = -X_0;
    end
end
 
function ret = sawtooth_wave(t,X_0, T_0)
    T = mod(t,T_0);     % T는 주기로 나눈 나머지
    ret = 2*X_0*(T/T_0) - X_0;
end
Colored by Color Scripter
cs

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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
clear;
% (ii) Angle plot를 출력
 
% function handle 선언
x_1 = @(t) square_wave(t, 1, 20);
x_2 = @(t) sawtooth_wave(t, 1, 0.02);
 
k = -10:1:10;
 
% 1. input : square wave
C_1 = my_fourier_coefficient(x_1, 20, k);       % Ck를 구한다.(k = -10:1:10)
t = 0:100/10000:100;
x_1_1 = zeros(1,10001);                         % 10001개의 x_1_1배열을 생성
for a = 1:1:10001
    x_1_1(a) = x_1(t(a));                       % x_1_1에 squere wave를 담아준다.
end
 
 
% 그래프 출력
figure(1); subplot(211);plot(t, x_1_1);
title('Square wave')
xlabel('t');
ylabel('x(t)');
figure(1); subplot(212); stem(k, angle(C_1)*180/pi);
title('Angle')
xlabel('k');
ylabel('degree of C_k');
 
 
 
% 2. input : sawtooth wave
C_2 = my_fourier_coefficient(x_2, 0.02, k);     % Ck를 구한다.(k = -10:1:10)
t = 0:0.1/10000:0.1;
x_2_2 = x_2(t);                                 % x_2_2에 sawtooth wave를 담아준다.
 
% 그래프 출력
figure(2); subplot(211);plot(t, x_2_2);
title('Sawtooth wave')
xlabel('t');
ylabel('x(t)');
figure(2); subplot(212); stem(k, angle(C_2)*180/pi);
title('Angle')
xlabel('k');
ylabel('degree of C_k');
 
 
 
function a = myintegral(fun, gap, min, max)
    % function handle 검사
    if ~isa(fun,'function_handle')
        error(message('MATLAB:integral:funArgNotHandle'));  % 오류메세지 출력 후 프로세스 종료
    end
    
    a = 0;                      % 결과로 반환될 a를 0으로 초기화
    
    for i = min : gap : max
        a = a + fun(i)*gap;     % 반복문 활용한 수치 적분
    end
end
 
% 결과를 배열에 담아 return 하는 convolution 함수
% x,h를 배열로 입력
function ret = my_fourier_coefficient(fun, T_0, k)
    % function handle 검사
    if ~isa(fun,'function_handle')
        error(message('MATLAB:integral:funArgNotHandle'));  % 오류메세지 출력 후 프로세스 종료
    end
    
 
    FUN = @(t) fun(t).*exp(-i*k*(2*pi/T_0).*t);             % 퓨리에 변환을 활용한다.
    ret = (1/T_0)*myintegral(FUN, T_0 / 10000, 0, T_0);     % 직접 만든 적분을 활용해 C_k을 구한다.
        
end
 
function ret = my_fourier_series(fun, N, T_0, t)
    % function handle 검사
    if ~isa(fun,'function_handle')
        error(message('MATLAB:integral:funArgNotHandle'));  % 오류메세지 출력 후 프로세스 종료
    end
    
    ret = 0;
    
    for k=-N:1:N
        C_k = my_fourier_coefficient(fun, T_0, k);
        ret = ret + C_k*exp(i*k*(2*pi/T_0)*t);
    end
end
 
function ret = square_wave(t,X_0, T_0)
    T = mod(t,T_0);     % T는 주기로 나눈 나머지
    if(T < (T_0 / 2))
        ret = X_0;
    else
        ret = -X_0;
    end
end
 
function ret = sawtooth_wave(t,X_0, T_0)
    T = mod(t,T_0);     % T는 주기로 나눈 나머지
    ret = 2*X_0*(T/T_0) - X_0;
end
Colored by Color Scripter
cs

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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
clear;
% (iii)
 
% function handle 선언
x_1 = @(t) square_wave(t, 1, 20);
x_2 = @(t) sawtooth_wave(t, 1, 0.02);
 
k = -10:1:10;
 
% 1. input : square wave
C_1 = my_fourier_coefficient(x_1, 20, k); 
t = 0:100/10000:100;
x_1_1 = zeros(1,10001);                 
for a = 1:1:10001
    x_1_1(a) = x_1(t(a));       
end
 
figure(1); subplot(311);plot(t, x_1_1);
title('Square wave')
xlabel('t');
ylabel('x(t)');
 
re_x_1_100 = my_fourier_series(x_1, 100, 20, t);
figure(1);subplot(312);plot(t, re_x_1_100);
title('Square wave in N = 100')
xlabel('t');
ylabel('x(t)');
 
re_x_1_10 = my_fourier_series(x_1, 10, 20, t);
figure(1);subplot(313);plot(t, re_x_1_10);
title('Square wave in N = 10')
xlabel('t');
ylabel('x(t)');
 
 
% 2. input : sawtooth wave
C_2 = my_fourier_coefficient(x_2, 0.02, k);     % Ck를 구한다.(k = -10:1:10)
t = 0:0.1/10000:0.1;
 
x_2_2 = x_2(t);                                 % x_2_2에 sawtooth wave를 담아준다.
figure(2); subplot(311);plot(t, x_2_2);
title('Sawtooth wave')
xlabel('t');
ylabel('x(t)');
 
re_x_2_100 = my_fourier_series(x_2, 100, 0.02, t);
figure(2);subplot(312);plot(t, re_x_2_100);
title('Sawtooth wave in N = 100')
xlabel('t');
ylabel('x(t)');
 
re_x_2_10 = my_fourier_series(x_2, 10, 0.02, t);
figure(2);subplot(313);plot(t, re_x_2_10);
title('Sawtooth wave in N = 10')
xlabel('t');
ylabel('x(t)');
 
 
 
 
function a = myintegral(fun, gap, min, max)
    % function handle 검사
    if ~isa(fun,'function_handle')
        error(message('MATLAB:integral:funArgNotHandle'));  % 오류메세지 출력 후 프로세스 종료
    end
    
    a = 0;                      % 결과로 반환될 a를 0으로 초기화
    
    for i = min : gap : max
        a = a + fun(i)*gap;     % 반복문 활용한 수치 적분
    end
end
 
% 결과를 배열에 담아 return 하는 convolution 함수
% x,h를 배열로 입력
function ret = my_fourier_coefficient(fun, T_0, k)
    % function handle 검사
    if ~isa(fun,'function_handle')
        error(message('MATLAB:integral:funArgNotHandle'));  % 오류메세지 출력 후 프로세스 종료
    end
    
 
    FUN = @(t) fun(t).*exp(-i*k*(2*pi/T_0).*t);             % 퓨리에 변환을 활용한다.
    ret = (1/T_0)*myintegral(FUN, T_0 / 10000, 0, T_0);     % 직접 만든 적분을 활용해 C_k을 구한다.
        
end
 
function ret = my_fourier_series(fun, N, T_0, t)
    % function handle 검사
    if ~isa(fun,'function_handle')
        error(message('MATLAB:integral:funArgNotHandle'));  % 오류메세지 출력 후 프로세스 종료
    end
    
    ret = 0;
    
    for k=-N:1:N
        C_k = my_fourier_coefficient(fun, T_0, k);
        ret = ret + C_k*exp(i*k*(2*pi/T_0)*t);
    end
end
 
function ret = square_wave(t,X_0, T_0)
    T = mod(t,T_0);     % T는 주기로 나눈 나머지
    if(T < (T_0 / 2))
        ret = X_0;
    else
        ret = -X_0;
    end
end
 
function ret = sawtooth_wave(t,X_0, T_0)
    T = mod(t,T_0);     % T는 주기로 나눈 나머지
    ret = 2*X_0*(T/T_0) - X_0;
end
Colored by Color Scripter
cs

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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
clear;
format long;
% (iv)
 
% function handle 선언
x_1 = @(t) square_wave(t, 1, 20);
x_2 = @(t) sawtooth_wave(t, 1, 0.02);
 
% 1. input : square wave
fun = @(t) abs(x_1(t))^2;
sol_1 = (1/(2*20))*myintegral(fun, 2*20 / 10000, -20, 20)
 
N = 10;
ret_1_N_10 = 0;
for a = -N:1:N
    ret_1_N_10 = ret_1_N_10 + abs(my_fourier_coefficient(x_1, 20, a))^2;
end
ret_1_N_10
 
 
N=100;
ret_1_N_100 = 0;
for a = -N:1:N
    ret_1_N_100 = ret_1_N_100 + abs(my_fourier_coefficient(x_1, 20, a))^2;
end
ret_1_N_100
 
% 2. input : sawtooth wave
fun = @(t) abs(x_2(t))^2;
sol_2 = (1/(2*0.02))*myintegral(fun, 2*0.02 / 10000, -0.02, 0.02)
 
N = 10;
ret_2_N_10 = 0;
for a = -N:1:N
    ret_2_N_10 = ret_2_N_10 + abs(my_fourier_coefficient(x_2, 0.02, a))^2;
end
ret_2_N_10
 
 
N=100;
ret_2_N_100 = 0;
for a = -N:1:N
    ret_2_N_100 = ret_2_N_100 + abs(my_fourier_coefficient(x_2, 0.02, a))^2;
end
ret_2_N_100
 
 
 
function a = myintegral(fun, gap, min, max)
    % function handle 검사
    if ~isa(fun,'function_handle')
        error(message('MATLAB:integral:funArgNotHandle'));  % 오류메세지 출력 후 프로세스 종료
    end
    
    a = 0;                      % 결과로 반환될 a를 0으로 초기화
    
    for i = min : gap : max
        a = a + fun(i)*gap;     % 반복문 활용한 수치 적분
    end
end
 
% 결과를 배열에 담아 return 하는 convolution 함수
% x,h를 배열로 입력
function ret = my_fourier_coefficient(fun, T_0, k)
    % function handle 검사
    if ~isa(fun,'function_handle')
        error(message('MATLAB:integral:funArgNotHandle'));  % 오류메세지 출력 후 프로세스 종료
    end
    
 
    FUN = @(t) fun(t).*exp(-i*k*(2*pi/T_0).*t);             % 퓨리에 변환을 활용한다.
    ret = (1/T_0)*myintegral(FUN, T_0 / 10000, 0, T_0);     % 직접 만든 적분을 활용해 C_k을 구한다.
        
end
 
function ret = my_fourier_series(fun, N, T_0, t)
    % function handle 검사
    if ~isa(fun,'function_handle')
        error(message('MATLAB:integral:funArgNotHandle'));  % 오류메세지 출력 후 프로세스 종료
    end
    
    ret = 0;
    
    for k=-N:1:N
        C_k = my_fourier_coefficient(fun, T_0, k);
        ret = ret + C_k*exp(i*k*(2*pi/T_0)*t);
    end
end
 
function ret = square_wave(t,X_0, T_0)
    T = mod(t,T_0);     % T는 주기로 나눈 나머지
    if(T < (T_0 / 2))
        ret = X_0;
    else
        ret = -X_0;
    end
end
 
function ret = sawtooth_wave(t,X_0, T_0)
    T = mod(t,T_0);     % T는 주기로 나눈 나머지
    ret = 2*X_0*(T/T_0) - X_0;
end
Colored by Color Scripter
cs

3. 참고

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

도움 되셨으면 하트 꾹!

저작자표시

'---------개인공부-------- > |matlab| Signals and Systems' 카테고리의 다른 글

[matlab] Discrete Fourier Transformation (DFT) (이산 퓨리에 변환!!!, 필터, 신호 및 시스템, 익명 함수, for, Signals and Systems) (0)	2020.07.19
[matlab] Continuous-Time Fourier Transformation (퓨리에 변환, 연속 시간, 적분, 신호 및 시스템, 익명 함수, Signals and Systems) (2)	2020.07.19
[matlab] Discrete Time Convolution (이산시간 컨볼루션, 신호 및 시스템, Signals and Systems) (0)	2020.07.19
[matlab] Plotting and Even and Odd Signals (그래프 출력, 신호 및 시스템, 홀함수 짝함수, Signals and Systems) (0)	2020.07.19
[matlab] Numerical Integration (적분, 신호 및 시스템, 익명 함수, for, Signals and Systems) (1)	2020.07.19

코린이의 작업공간

[matlab] Continuous-time Fourier Series (연속 시간, 퓨리에 시리즈, function handle 활용, 함수의 신기한 활용,신호 및 시스템, Signals and Systems)

'---------개인공부-------- > |matlab| Signals and Systems' 카테고리의 다른 글

+ Recent posts

티스토리툴바