-
Notifications
You must be signed in to change notification settings - Fork 1
/
my_conv2.m
72 lines (62 loc) · 1.8 KB
/
my_conv2.m
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
function h = my_conv2(f, g, conv_type)
% FFTタイプの畳み込み
% 多項式 f, g の積を計算する。
%
% Parameters
% ----------
% f : f[i] に、x^i の係数が入っている
%
% g : g[i] に、x^i の係数が入っている
%
%
% Returns
% -------
% h : f,g の積
%
% h の長さ以上の n=2^k を計算
arguments
f = []
g = []
conv_type = 'full'
end
%% サイズ調整
fft_len_x = 1;
fft_len_y = 1;
while (2 * fft_len_x) < (size(f,2) + size(g,2) - 1)
fft_len_x = fft_len_x * 2;
end
fft_len_x = fft_len_x * 2 ;
while (2 * fft_len_y) < (size(f,1) + size(g,1) - 1)
fft_len_y = fft_len_y * 2;
end
fft_len_y = fft_len_y * 2 ;
% フーリエ変換で実装する
Ff = fft2(f, fft_len_y,fft_len_x);% Ff = real(Ff);
Fg = fft2(g, fft_len_y,fft_len_x);% Fg = real(Fg);
% 各点積
Fh = Ff .* Fg;
% フーリエ逆変換
h = ifft2(Fh);
% フルサイズの取得
len_x = size(f,2) + size(g,2) - 1;
len_y = size(f,1) + size(g,1) - 1;
h = h(1:1+len_y-1, 1:1+len_x-1);
% same, valid, fullで分ける
switch conv_type
case 'same'
len_x = size(f,2);
len_y = size(f,1);
case 'valid'
len_x = max(size(f,2)-size(g,2) + 1, 0);
len_y = max(size(f,1)-size(g,1) + 1, 0);
% case 'full'
% len_x = size(f,2) + size(g,2) - 1;
% len_y = size(f,1) + size(g,1) - 1;
% otherwise
% len_x = size(f,2) + size(g,2) - 1;
% len_y = size(f,1) + size(g,1) - 1;
end
st_x = ceil((size(h, 2) - len_x) / 2) + 1;
st_y = ceil((size(h, 1) - len_y) / 2) + 1;
h = (h(st_y:st_y + len_y - 1, st_x:st_x + len_x - 1));
end