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:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
152:
153:
154:
155:
156:
157:
158:
159:
160:
161:
162:
163:
164:
165:
166:
167:
168:
169:
170:
171:
172:
173:
174:
175:
176:
177:
178:
179:
180:
181:
182:
183:
184:
185:
186:
187:
188:
189:
190:
191:
192:
193:
194:
195:
196:
197:
198:
199:
200:
201:
202:
203:
204:
205:
206:
207:
208:
209:
210:
211:
212:
213:
214:
215:
216:
217:
218:
219:
220:

%trigonometrische Interpolation
clear all;

disp(' ');
disp('Es wird trigonometrisch approximiert/interpoliert.');
func=input('Vergleichsfunktion in trigof anlegen? 0 -ja, 1 - nein: ');
disp(' ');
disp('Der Datensatz hat die Form');
disp(' Knoten: [x_0,...,x_n]');
disp(' f-Werte: [y_0,...,y_n]');
disp(' ');
%yi=input('f-Werte eingeben: ');
n=16;%length(yi);
nexp=log(n)/log(2);
x=linspace(0,2*pi);
disp('---------------------------------------------------');

%0. Einheitswurzeln erzeigen
for k=1:n
xi(1,k)=2*pi*(k-1)/n;
yi(1,k)=trigof(xi(1,k));
end
yi(1,1)=0.5*(2*pi)^0.5;



%1. Funktion f plotten
if func==0
%abschnittsweise Funktion f anlegen.
y=sqrt(x);

%Plot erstellen
subplot(4,1,1)
plot(xi,yi,'ro',x,y),
%plot(x,y,xi,yi,'ro'),
axis([0,2*pi,0,3]),
title('Funktion f'),
hold;

%2. Aproximation

disp(' ');
disp('1. Approximation. Koeffizienten a,b');
disp(' ');
%Integrale für die Koeffizienten werden mit summierter Trapezregel bestimmt (Genauigkeit eps)
%Berechnungen der Koeffizienten a und b
eps=0.01;
maxd=1;

N=ceil(sqrt(2*pi^3*(n-1)^2/(3*eps)));

%Knoten für summierte Trapezregel bestimmen
for m=1:N
xs(1,m)=2*pi/N*(m-1);
end

sum13=0;
for j=0:n
sum11=0;
sum12=0;

for k=1:N
sum11=sum11+trigof(xs(1,k)).*sin(j*xs(1,k));
sum12=sum12+trigof(xs(1,k)).*cos(j*xs(1,k));
end
a(j+1,1)=2/N*sum12;
b(j+1,1)=2/N*sum11;

if j>0
sum13=sum13+a(j+1,1).*cos(j*x)+b(j+1,1).*sin(j*x);
end
end
sum13=sum13+0.5*a(1,1);

TNA=sum13;
a
b

%Plot erstellen
subplot(4,1,2)
plot(x,TNA,xi,yi,'ro'),
axis([0,2*pi,0,3]),
title('Fourier-Polynom');
end



%3. Trigonomatrische Interpolation
disp(' ');
disp('2. Trigonometrische Interpolation. Koeffizienten d');
disp(' ')
%Integrale
%Berechnung der Einheitswurzeln
R=1; phi=(2*pi)/n;
w=R.*exp(i*phi);

%Berechnung der komplexen Koeffizienten d

if floor(nexp)==ceil(nexp)

%Mit FFT
yff=yi';
d=1/n.*fft(yff);
d=bitrevorder(d);

else

%Ohne FFT
for j=0:n-1
sum31=0;
for k=0:n-1
sum31=sum31+yi(1,k+1).*w^(-j*k);
end
d(j+1,1)=(1/n).*sum31;
end
d;
end
d

%Trigonometriches IP bestimmen und plotfähig machen
sum32=0;
for j=0:n-1
sum32=sum32 + d(j+1,1).*(cos(j*x)+i.*sin(j*x));
end
TN2=sum32;

%Realteil
TN2R=real(TN2);

%Imaginärteil
TN2C=imag(TN2);

%Plot erstellen
subplot(4,1,3)
plot(x,TN2R,x,TN2C,xi,yi,'ro'),
axis([0,2*pi,0,3]),
title('Trigonometrische Interpolation');
legend('Realteil','Imaginärteil'),
hold


%4. Diskrete Fourier-Analyse
disp(' ')
disp('3. Diskrete Fourier Analyse. Koeffizienten ai, bi');
disp(' ')
%Berechnung von m
if 0==mod(n,2)
m=n/2;
else
m=(n+1)/2;
end
m;

%Berechnungen der Koeffizienten a und b
if 0==mod(n,2)
for j=0:m
sum1=0;
for k=0:n-1
sum1=sum1+yi(1,k+1).*cos(k*xi(1,j+1));
end
ai(j+1,1)=2/n*sum1;
end
ai
else
for j=0:m-1
sum1=0;
for k=0:n-1
sum1=sum1+yi(1,k+1).*cos(k*xi(1,j+1));
end
ai(j+1,1)=2/n*sum1;
end
ai
end

for j=1:m-1
sum2=0;
for k=0:n-1
sum2=sum2+yi(1,k+1).*sin(k.*xi(1,j+1));
end
bi(j,1)=2/n*sum2;
end
bi

%DFA-Funktion bestimmen und plotfähig machen
sum3=0;
for j=1:m-1
sum3=sum3 + (ai(j+1,1).*cos(j*x)+bi(j,1).*sin(j*x));
end

if 0==mod(n,2)
TNDFA=ai(1,1)./2+sum3+ai(m+1,1)./2.*cos(m*x);
else
TNDFA=ai(1,1)./2+sum3;
end

%abschnittsweise Funktion f anlegen.
x1=0;
x2=linspace(0,pi);
x3=pi;
x4=linspace(pi,2*pi);

fx1=0;
for q=1:length(x2)
fx2(1,q)=-1;
end
fx3=0;
for j=1:length(x4)
fx4(1,q)=1;
end


%Plot erstellen
subplot(4,1,4)
plot(x,TNDFA,xi,yi,'ro'),
axis([0,2*pi,0,3]),
title('Diskrete Fourieranalyse');
hold