-
Notifications
You must be signed in to change notification settings - Fork 0
/
myfft32.m
129 lines (81 loc) · 2.18 KB
/
myfft32.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
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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
function X = myfft32(x)
%% calculate the DFT of a complex vector of length 32
%verify length
if length(x)~=32
'error-input length is not 32'
return
end
%% define W_N
N_16 = 16;
W_16 = exp((-j*2*pi)/N_16);
N_32=32;
W_32 = exp((-j*2*pi)/N_32);
%% separate odds from evens into TWO 16-SAMPLE arrays; Keep in mind matlab has it all backwards
y=zeros(1,16);
z=zeros(1,16);
% essentially picking every other value for each array
for n = 0:15
y(n+1)=x(2*n+1); % pick up EVEN subscripts
z(n+1)=x(2*n+2); % pick up ODD subscripts
end
%% separeate further into FOUR 8-SAMPLE arrays;
yy0 = zeros(1,8);
yy2 = zeros(1,8);
zz1 = zeros(1,8);
zz3 = zeros(1,8);
for n = 0:7
yy0(n+1)=y(2*n+1); % pick up EVEN subscripts
yy2(n+1)=y(2*n+2); % pick up ODD subscripts
end
for n = 0:7
zz1(n+1)=z(2*n+1); % pick up EVEN subscripts
zz3(n+1)=z(2*n+2); % pick up ODD subscripts
end
%% calculaing EVEN-coefficient FFTs
YY0 = fft(yy0); % FFT of the sequence starting with subscript 0
YY2 = fft(yy2); % FFT of the sequence starting with subscript 2
Y = zeros(1,16);
for n=1:8
Y(n) = YY0(n)+ W_16^(n-1)*YY2(n);
end
for n = 9:16
Y(n) = YY0(n-8)- W_16^(n-9)*YY2(n-8);
end
%% calculaing ODD-coefficient FFTs
ZZ1 = fft(zz1); % FFT of the sequence starting with subscript 1
ZZ3 = fft(zz3); % FFT of the sequence starting with subscript 3
Z = zeros(1,16);
for n=1:8
Z(n) = ZZ1(n)+ W_16^(n-1)*ZZ3(n);
end
for n = 9:16
Z(n) = ZZ1(n-8)- W_16^(n-9)*ZZ3(n-8);
end
%% put it all together
X=zeros(1,32);
for n=1:16
X(n) = Y(n)+ W_32^(n-1)*Z(n);
end
for n = 17:32
X(n) = Y(n-16)- W_32^(n-17)*Z(n-16);
end
%% Normalize the calculated DFT by dividing by N
X=X/length(X);
%% Plot
figure()
stem(fftshift(abs(X)),'r');
%legend('My FFT implementation');
hold on
stem(fftshift(abs(fft(x)/length(fft(x)))),'b');
%legend('MATLAB FFT reference');
hold off
legend('My FFT implementation','MATLAB FFT reference');
title("fftshift(abs(X/N)) & fftshift(abs(fft(x)/N))")
figure()
subplot(2,1,1)
stem(fftshift(abs(X)),'r');
legend('My FFT implementation');
subplot(2,1,2)
stem(fftshift(abs(fft(x)/length(fft(x)))),'b');
legend('MATLAB FFT reference');
end