|
Python |
Matlab/Octave |
| Vorbereitung (Laden von Modulen/Paketen) |
#python -m pip install --upgrade numpy
#python -m pip install --upgrade matplotlib
#pip install --upgrade numpy
#pip install --upgrade matplotlib
from numpy import *
from numpy.random import *
from matplotlib.pyplot import *
|
%Matlab: erfordert die Statistics and Machine Learning Toolbox
%Octave: pkg load statistics
if exist('OCTAVE_VERSION','builtin')~=0
pkg load statistics
end
|
| Signalrekonstruktion (Ausschnitt eines AR1-Prozess) |
dt=1.0
N=10
T=N*dt
Np=50
phi=0.9
mx=0
vx=1
theta=sqrt((1-phi**2)*vx)
tp=arange(0,Np)*dt-2*T
xp=zeros(Np)
xp[0]=normal(mx,sqrt(vx))
for i in range(1,Np):
xp[i]=phi*(xp[i-1]-mx)+normal(mx,theta)
t=arange(0,N)*dt
x=xp[(Np-N)//2:(Np-N)//2+N]
dtrek=0.1
trek=arange(0,int(round(3*T/dtrek)))*dtrek-T
xrek=zeros(len(trek))
for j in range(0,len(xrek)):
for i in range(0,len(xp)):
xrek[j]+=xp[i]*sinc(trek[j]/float(dt)-tp[i])
plot(t,x,'o',trek,xrek)
show()
|
dt=1.0;
N=10;
T=N*dt;
Np=50;
phi=0.9;
mx=0;
vx=1;
theta=sqrt((1-phi^2)*vx);
tp=(0:Np-1)*dt-2*T;
xp=zeros(1,Np);
xp(1)=normrnd(mx,sqrt(vx));
for i=2:Np
xp(i)=phi*(xp(i-1)-mx)+normrnd(mx,theta);
end
t=(0:N-1)*dt;
x=xp(floor((Np-N)/2)+1:floor((Np-N)/2)+N);
dtrek=0.1;
trek=(0:fix(round(3*T/dtrek))-1)*dtrek-T;
xrek=zeros(1,length(trek));
for j=1:length(xrek)
for i=1:length(xp)
xrek(j)=xrek(j)+xp(i)*sinc(trek(j)/dt-tp(i));
end
end
plot(t,x,'o',trek,xrek)
|
| Signalrekonstruktion (Interpretation als Periode eines periodischen Signals und zeitlich begrenztes Energiesignal) |
dtrek=0.1
trek=arange(0,int(round(3*T/dtrek)))*dtrek-T
xrek=zeros(len(trek))
for j in range(0,len(xrek)):
if abs(trek[j]/dt-round(trek[j]/dt))>1e-9:
if len(x)%2==0:
for i in range(0,len(x)):
xrek[j]+=x[i]*sin(pi*(trek[j]-t[i])/dt)*cos(pi*(trek[j]-t[i])/dt/float(len(x)))/sin(pi*(trek[j]-t[i])/dt/float(len(x)))/float(len(x))
else:
for i in range(0,len(x)):
xrek[j]+=x[i]*sin(pi*(trek[j]-t[i])/dt)/sin(pi*(trek[j]-t[i])/dt/float(len(x)))/float(len(x))
else:
xrek[j]=x[int(round(trek[j]/dt))%len(x)]
plot(t,x,'o',trek,xrek)
show()
|
dtrek=0.1;
trek=(0:fix(round(3*T/dtrek))-1)*dtrek-T;
xrek=zeros(1,length(trek));
for j=1:length(xrek)
if abs(trek(j)/dt-fix(trek(j)/dt))>1e-9
if mod(length(x),2)==0
for i=1:length(x)
xrek(j)=xrek(j)+x(i)*sin(pi*(trek(j)-t(i))/dt)*cos(pi*(trek(j)-t(i))/dt/length(x))/sin(pi*(trek(j)-t(i))/dt/length(x))/length(x);
end
else
for j=1:length(x)
xrek(j)=xrek(j)+x(i)*sin(pi*(trek(j)-t(i))/dt)/sin(pi*(trek(j)-t(i))/dt/length(x))/length(x);
end
end
else
xrek(j)=x(mod(fix(trek(j)/dt),length(x))+1);
end
end
plot(t,x,'o',trek,xrek)
|
| Signalrekonstruktion (Interpretation als zeitlich begrenztes Energiesignal) |
dtrek=0.1
trek=arange(0,int(round(3*T/dtrek)))*dtrek-T
xrek=zeros(len(trek))
for j in range(0,len(xrek)):
for i in range(0,len(x)):
xrek[j]+=x[i]*sinc(trek[j]/float(dt)-t[i])
plot(t,x,'o',trek,xrek,append(arange(-int(round(T/float(dt))),0),arange(int(round(T/float(dt))),int(round(2*T/float(dt))))),zeros(int(round(2*T/float(dt)))),'.')
show()
|
dtrek=0.1;
trek=(0:fix(round(3*T/dtrek))-1)*dtrek-T;
xrek=zeros(1,length(trek));
for j=1:length(xrek)
for i=1:length(x)
xrek(j)=xrek(j)+x(i)*sinc(trek(j)/dt-t(i));
end
end
plot(t,x,'o',trek,xrek)
|
| bandbegrenztes Signal, Abtastung unter Einhaltung des Abtasttheorems |
T=100.0
f=1/20.0
dtsuper=0.1
tsuper=arange(0,int(round(T/dtsuper)))*dtsuper
xsuper=sin(2*pi*tsuper*f)
dtabt=1.0
tabt=arange(0,int(round(T/dtabt)))*dtabt
xabt=sin(2*pi*tabt*f)
plot(tsuper,xsuper,tabt,xabt,'o')
show()
|
T=100;
f=1/20;
dtsuper=0.1;
tsuper=(0:fix(round(T/dtsuper))-1)*dtsuper;
xsuper=sin(2*pi*tsuper*f);
dtabt=1;
tabt=(0:fix(round(T/dtabt))-1)*dtabt;
xabt=sin(2*pi*tabt*f);
plot(tsuper,xsuper,tabt,xabt,'o')
|
| Leistungsspektrum |
from numpy.fft import *
N=int(round(T/dtsuper))
fsuper=roll(arange(-(N//2),(N+1)//2)/float(N*dtsuper),(N+1)//2)
X=fft(xsuper)
Psuper=abs(X)**2/N**2
N=int(round(T/dtabt))
fabt=roll(arange(-(N//2),(N+1)//2)/float(N*dtabt),(N+1)//2)
X=fft(xabt)
Pabt=abs(X)**2/N**2
vlines(fsuper,0,Psuper)
plot(fabt,Pabt,'.')
xlim(-0.25,0.25)
show()
|
N=fix(round(T/dtsuper));
fsuper=(mod((0:N-1)+fix(N/2),N)-fix(N/2))/(N*dtsuper);
X=fft(xsuper);
Psuper=abs(X).^2/N^2;
N=fix(round(T/dtabt));
fabt=(mod((0:N-1)+fix(N/2),N)-fix(N/2))/(N*dtabt);
X=fft(xabt);
Pabt=abs(X).^2/N^2;
stem(fsuper,Psuper,'.')
hold on
plot(fabt,Pabt,'o')
hold off
xlim([-0.25 0.25])
|
| vollständige Rekonstruktion |
dtrek=0.1
trek=arange(0,int(round(T/dtrek)))*dtrek
xrek=zeros(len(trek))
for j in range(0,len(xrek)):
if abs(trek[j]/dtabt-round(trek[j]/dtabt))>1e-9:
if len(xabt)%2==0:
for i in range(0,len(xabt)):
xrek[j]+=xabt[i]*sin(pi*(trek[j]-tabt[i])/dtabt)*cos(pi*(trek[j]-tabt[i])/dtabt/float(len(xabt)))/sin(pi*(trek[j]-tabt[i])/dtabt/float(len(xabt)))/float(len(xabt))
else:
for i in range(0,len(xabt)):
xrek[j]+=xabt[i]*sin(pi*(trek[j]-tabt[i])/dtabt)/sin(pi*(trek[j]-tabt[i])/dtabt/float(len(xabt)))/float(len(xabt))
else:
xrek[j]=xabt[int(round(trek[j]/dtabt))%len(xabt)]
plot(tsuper,xsuper,tabt,xabt,'o')
plot(trek,xrek,'--',dashes=[10,10])
show()
|
dtrek=0.1;
trek=(0:fix(round(T/dtrek))-1)*dtrek;
xrek=zeros(1,length(trek));
for j=1:length(xrek)
if abs(trek(j)/dtabt-fix(trek(j)/dtabt))>1e-9
if mod(length(xabt),2)==0
for i=1:length(xabt)
xrek(j)=xrek(j)+xabt(i)*sin(pi*(trek(j)-tabt(i))/dtabt)*cos(pi*(trek(j)-tabt(i))/dtabt/length(xabt))/sin(pi*(trek(j)-tabt(i))/dtabt/length(xabt))/length(xabt);
end
else
for i=1:length(xabt)
xrek(j)=xrek(j)+xabt(i)*sin(pi*(trek(j)-tabt(i))/dtabt)/sin(pi*(trek(j)-tabt(i))/dtabt/length(xabt))/length(xabt);
end
end
else
xrek(j)=xabt(mod(fix(trek(j)/dtabt),length(xabt))+1);
end
end
plot(tsuper,xsuper,tabt,xabt,'o',trek,xrek,'--')
|
| Unterabtastung |
T=100.0
f=17.0/100.0
dtsuper=0.1
tsuper=arange(0,int(round(T/dtsuper)))*dtsuper
xsuper=sin(2*pi*tsuper*f)
dtabt=5.0
tabt=arange(0,int(round(T/dtabt)))*dtabt
xabt=sin(2*pi*tabt*f)
plot(tsuper,xsuper,tabt,xabt,'o')
show()
|
T=100;
f=17/100;
dtsuper=0.1;
tsuper=(0:fix(round(T/dtsuper))-1)*dtsuper;
xsuper=sin(2*pi*tsuper*f);
dtabt=5;
tabt=(0:fix(round(T/dtabt))-1)*dtabt;
xabt=sin(2*pi*tabt*f);
plot(tsuper,xsuper,tabt,xabt,'o')
|
| Leistungsspektrum |
from numpy.fft import *
N=int(round(T/dtsuper))
fsuper=roll(arange(-(N//2),(N+1)//2)/float(N*dtsuper),(N+1)//2)
X=fft(xsuper)
Psuper=abs(X)**2/N**2
N=int(round(T/dtabt))
fabt=roll(arange(-(N//2),(N+1)//2)/float(N*dtabt),(N+1)//2)
X=fft(xabt)
Pabt=abs(X)**2/N**2
vlines(fsuper,0,Psuper)
plot(fabt,Pabt,'.')
xlim(-0.25,0.25)
show()
|
N=fix(round(T/dtsuper));
fsuper=(mod((0:N-1)+fix(N/2),N)-fix(N/2))/(N*dtsuper);
X=fft(xsuper);
Psuper=abs(X).^2/N^2;
N=fix(round(T/dtabt));
fabt=(mod((0:N-1)+fix(N/2),N)-fix(N/2))/(N*dtabt);
X=fft(xabt);
Pabt=abs(X).^2/N^2;
stem(fsuper,Psuper,'.')
hold on
plot(fabt,Pabt,'o')
hold off
xlim([-0.25 0.25])
|
| fehlerhafte Rekonstruktion trotz korrekter Abtastpunkte |
dtrek=0.1
trek=arange(0,int(round(T/dtrek)))*dtrek
xrek=zeros(len(trek))
for j in range(0,len(xrek)):
if abs(trek[j]/dtabt-round(trek[j]/dtabt))>1e-9:
if len(xabt)%2==0:
for i in range(0,len(xabt)):
xrek[j]+=xabt[i]*sin(pi*(trek[j]-tabt[i])/dtabt)*cos(pi*(trek[j]-tabt[i])/dtabt/float(len(xabt)))/sin(pi*(trek[j]-tabt[i])/dtabt/float(len(xabt)))/float(len(xabt))
else:
for i in range(0,len(xabt)):
xrek[j]+=xabt[i]*sin(pi*(trek[j]-tabt[i])/dtabt)/sin(pi*(trek[j]-tabt[i])/dtabt/float(len(xabt)))/float(len(xabt))
else:
xrek[j]=xabt[int(round(trek[j]/dtabt))%len(xabt)]
plot(tsuper,xsuper,tabt,xabt,'o')
plot(trek,xrek,'--',dashes=[10,10])
show()
|
dtrek=0.1;
trek=(0:fix(round(T/dtrek))-1)*dtrek;
xrek=zeros(1,length(trek));
for j=1:length(xrek)
if abs(trek(j)/dtabt-fix(trek(j)/dtabt))>1e-9
if mod(length(xabt),2)==0
for i=1:length(xabt)
xrek(j)=xrek(j)+xabt(i)*sin(pi*(trek(j)-tabt(i))/dtabt)*cos(pi*(trek(j)-tabt(i))/dtabt/length(xabt))/sin(pi*(trek(j)-tabt(i))/dtabt/length(xabt))/length(xabt);
end
else
for i=1:length(xabt)
xrek(j)=xrek(j)+xabt(i)*sin(pi*(trek(j)-tabt(i))/dtabt)/sin(pi*(trek(j)-tabt(i))/dtabt/length(xabt))/length(xabt);
end
end
else
xrek(j)=xabt(mod(fix(trek(j)/dtabt),length(xabt))+1);
end
end
plot(tsuper,xsuper,tabt,xabt,'o',trek,xrek,'--')
|
| Unterabtastung in Kombination mit nicht vollständig erfassten Perioden |
T=100.0
f=1/8.0
dtsuper=0.1
tsuper=arange(0,int(round(T/dtsuper)))*dtsuper
xsuper=sin(2*pi*tsuper*f)
dtabt=5.0
tabt=arange(0,int(round(T/dtabt)))*dtabt
xabt=sin(2*pi*tabt*f)
plot(tsuper,xsuper,tabt,xabt,'o')
show()
|
T=100;
f=1/8;
dtsuper=0.1;
tsuper=(0:fix(round(T/dtsuper))-1)*dtsuper;
xsuper=sin(2*pi*tsuper*f);
dtabt=5;
tabt=(0:fix(round(T/dtabt))-1)*dtabt;
xabt=sin(2*pi*tabt*f);
plot(tsuper,xsuper,tabt,xabt,'o')
|
| Leistungsspektrum |
from numpy.fft import *
N=int(round(T/dtsuper))
fsuper=roll(arange(-(N//2),(N+1)//2)/float(N*dtsuper),(N+1)//2)
X=fft(xsuper)
Psuper=abs(X)**2/N**2
N=int(round(T/dtabt))
fabt=roll(arange(-(N//2),(N+1)//2)/float(N*dtabt),(N+1)//2)
X=fft(xabt)
Pabt=abs(X)**2/N**2
vlines(fsuper,0,Psuper)
plot(fabt,Pabt,'.')
xlim(-0.25,0.25)
show()
|
N=fix(round(T/dtsuper));
fsuper=(mod((0:N-1)+fix(N/2),N)-fix(N/2))/(N*dtsuper);
X=fft(xsuper);
Psuper=abs(X).^2/N^2;
N=fix(round(T/dtabt));
fabt=(mod((0:N-1)+fix(N/2),N)-fix(N/2))/(N*dtabt);
X=fft(xabt);
Pabt=abs(X).^2/N^2;
stem(fsuper,Psuper,'.')
hold on
plot(fabt,Pabt,'o')
hold off
xlim([-0.25 0.25])
|
| fehlerhafte Rekonstruktion trotz korrekter Abtastpunkte |
dtrek=0.1
trek=arange(0,int(round(T/dtrek)))*dtrek
xrek=zeros(len(trek))
for j in range(0,len(xrek)):
if abs(trek[j]/dtabt-round(trek[j]/dtabt))>1e-9:
if len(xabt)%2==0:
for i in range(0,len(xabt)):
xrek[j]+=xabt[i]*sin(pi*(trek[j]-tabt[i])/dtabt)*cos(pi*(trek[j]-tabt[i])/dtabt/float(len(xabt)))/sin(pi*(trek[j]-tabt[i])/dtabt/float(len(xabt)))/float(len(xabt))
else:
for i in range(0,len(xabt)):
xrek[j]+=xabt[i]*sin(pi*(trek[j]-tabt[i])/dtabt)/sin(pi*(trek[j]-tabt[i])/dtabt/float(len(xabt)))/float(len(xabt))
else:
xrek[j]=xabt[int(round(trek[j]/dtabt))%len(xabt)]
plot(tsuper,xsuper,tabt,xabt,'o')
plot(trek,xrek,'--',dashes=[10,10])
show()
|
dtrek=0.1;
trek=(0:fix(round(T/dtrek))-1)*dtrek;
xrek=zeros(1,length(trek));
for j=1:length(xrek)
if abs(trek(j)/dtabt-fix(trek(j)/dtabt))>1e-9
if mod(length(xabt),2)==0
for i=1:length(xabt)
xrek(j)=xrek(j)+xabt(i)*sin(pi*(trek(j)-tabt(i))/dtabt)*cos(pi*(trek(j)-tabt(i))/dtabt/length(xabt))/sin(pi*(trek(j)-tabt(i))/dtabt/length(xabt))/length(xabt);
end
else
for i=1:length(xabt)
xrek(j)=xrek(j)+xabt(i)*sin(pi*(trek(j)-tabt(i))/dtabt)/sin(pi*(trek(j)-tabt(i))/dtabt/length(xabt))/length(xabt);
end
end
else
xrek(j)=xabt(mod(fix(trek(j)/dtabt),length(xabt))+1);
end
end
plot(tsuper,xsuper,tabt,xabt,'o',trek,xrek,'--')
|
| Frequenzgemisch |
T=100.0
f1=1/20.9
f2=1/6.85
dtsuper=0.1
tsuper=arange(0,int(round(T/dtsuper)))*dtsuper
xsuper=sin(2*pi*tsuper*f1+pi/3)+sin(2*pi*tsuper*f2)
dtabt=5.0
tabt=arange(0,int(round(T/dtabt)))*dtabt
xabt=sin(2*pi*tabt*f1+pi/3)+sin(2*pi*tabt*f2)
plot(tsuper,xsuper,tabt,xabt,'o')
show()
|
T=100;
f1=1/20.9;
f2=1/6.85;
dtsuper=0.1;
tsuper=(0:fix(round(T/dtsuper))-1)*dtsuper;
xsuper=sin(2*pi*tsuper*f1+pi/3)+sin(2*pi*tsuper*f2);
dtabt=5;
tabt=(0:fix(round(T/dtabt))-1)*dtabt;
xabt=sin(2*pi*tabt*f1+pi/3)+sin(2*pi*tabt*f2);
plot(tsuper,xsuper,tabt,xabt,'o')
|
| Leistungsspektrum |
from numpy.fft import *
N=int(round(T/dtsuper))
fsuper=roll(arange(-(N//2),(N+1)//2)/float(N*dtsuper),(N+1)//2)
X=fft(xsuper)
Psuper=abs(X)**2/N**2
N=int(round(T/dtabt))
fabt=roll(arange(-(N//2),(N+1)//2)/float(N*dtabt),(N+1)//2)
X=fft(xabt)
Pabt=abs(X)**2/N**2
vlines(fsuper,0,Psuper)
plot(fabt,Pabt,'.')
xlim(-0.25,0.25)
show()
|
N=fix(round(T/dtsuper));
fsuper=(mod((0:N-1)+fix(N/2),N)-fix(N/2))/(N*dtsuper);
X=fft(xsuper);
Psuper=abs(X).^2/N^2;
N=fix(round(T/dtabt));
fabt=(mod((0:N-1)+fix(N/2),N)-fix(N/2))/(N*dtabt);
X=fft(xabt);
Pabt=abs(X).^2/N^2;
stem(fsuper,Psuper,'.')
hold on
plot(fabt,Pabt,'o')
hold off
xlim([-0.25 0.25])
|
| fehlerhafte Rekonstruktion beider Frequenzanteile trotz korrekter Abtastpunkte |
dtrek=0.1
trek=arange(0,int(round(T/dtrek)))*dtrek
xrek=zeros(len(trek))
for j in range(0,len(xrek)):
if abs(trek[j]/dtabt-round(trek[j]/dtabt))>1e-9:
if len(xabt)%2==0:
for i in range(0,len(xabt)):
xrek[j]+=xabt[i]*sin(pi*(trek[j]-tabt[i])/dtabt)*cos(pi*(trek[j]-tabt[i])/dtabt/float(len(xabt)))/sin(pi*(trek[j]-tabt[i])/dtabt/float(len(xabt)))/float(len(xabt))
else:
for i in range(0,len(xabt)):
xrek[j]+=xabt[i]*sin(pi*(trek[j]-tabt[i])/dtabt)/sin(pi*(trek[j]-tabt[i])/dtabt/float(len(xabt)))/float(len(xabt))
else:
xrek[j]=xabt[int(round(trek[j]/dtabt))%len(xabt)]
plot(tsuper,xsuper,tabt,xabt,'o')
plot(trek,xrek,'--',dashes=[10,10])
show()
|
dtrek=0.1;
trek=(0:fix(round(T/dtrek))-1)*dtrek;
xrek=zeros(1,length(trek));
for j=1:length(xrek)
if abs(trek(j)/dtabt-fix(trek(j)/dtabt))>1e-9
if mod(length(xabt),2)==0
for i=1:length(xabt)
xrek(j)=xrek(j)+xabt(i)*sin(pi*(trek(j)-tabt(i))/dtabt)*cos(pi*(trek(j)-tabt(i))/dtabt/length(xabt))/sin(pi*(trek(j)-tabt(i))/dtabt/length(xabt))/length(xabt);
end
else
for i=1:length(xabt)
xrek(j)=xrek(j)+xabt(i)*sin(pi*(trek(j)-tabt(i))/dtabt)/sin(pi*(trek(j)-tabt(i))/dtabt/length(xabt))/length(xabt);
end
end
else
xrek(j)=xabt(mod(fix(trek(j)/dtabt),length(xabt))+1);
end
end
plot(tsuper,xsuper,tabt,xabt,'o',trek,xrek,'--')
|
| Anti-Aliasing-Filter vor der Abtastung |
from numpy.fft import *
T=100.0
f1=1/20.9
f2=1/6.85
dtsuper=0.1
tsuper=arange(0,int(round(T/dtsuper)))*dtsuper
xsuper=sin(2*pi*tsuper*f1)+0.5*sin(2*pi*tsuper*f2)
N=int(round(T/dtsuper))
fsuper=roll(arange(-(N//2),(N+1)//2)/float(N*dtsuper),(N+1)//2)
X=fft(xsuper)
dtabt=5.0
for i in range(0,len(X)):
if abs(fsuper[i])>=1/float(2*dtabt):
X[i]=0.0
xfilt=real(ifft(X))
tabt=arange(0,int(round(T/dtabt)))*dtabt
xabt=zeros(len(tabt))
for i in range(0,len(xabt)):
for j in range(0,len(xsuper)):
xabt[i]+=xfilt[j]*(sinc((tabt[i]-tsuper[j])/dtsuper)+sinc((tabt[i]-tsuper[j]+T)/dtsuper)+sinc((tabt[i]-tsuper[j]-T)/dtsuper))
plot(tsuper,xsuper,tsuper,xfilt,tabt,xabt,'o')
show()
|
T=100;
f1=1/20.9;
f2=1/6.85;
dtsuper=0.1;
tsuper=(0:fix(round(T/dtsuper))-1)*dtsuper;
xsuper=sin(2*pi*tsuper*f1+pi/3)+sin(2*pi*tsuper*f2);
N=fix(round(T/dtsuper));
fsuper=(mod((0:N-1)+fix(N/2),N)-fix(N/2))/(N*dtsuper);
X=fft(xsuper);
dtabt=5;
for i=1:length(X)
if abs(fsuper(i))>=1/(2*dtabt)
X(i)=0;
end
end
xfilt=real(ifft(X));
tabt=(0:fix(round(T/dtabt))-1)*dtabt;
xabt=zeros(1,length(tabt));
for i=1:length(xabt)
for j=1:length(xsuper)
xabt(i)=xabt(i)+xfilt(j)*(sinc((tabt(i)-tsuper(j))/dtsuper)+sinc((tabt(i)-tsuper(j)+T)/dtsuper)+sinc((tabt(i)-tsuper(j)-T)/dtsuper));
end
end
plot(tsuper,xsuper,tsuper,xfilt,tabt,xabt,'o')
|
| Leistungsspektrum |
from numpy.fft import *
N=int(round(T/dtsuper))
fsuper=roll(arange(-(N//2),(N+1)//2)/float(N*dtsuper),(N+1)//2)
X=fft(xfilt)
Pfilt=abs(X)**2/N**2
N=int(round(T/dtabt))
fabt=roll(arange(-(N//2),(N+1)//2)/float(N*dtabt),(N+1)//2)
X=fft(xabt)
Pabt=abs(X)**2/N**2
vlines(fsuper,0,Pfilt)
plot(fabt,Pabt,'.')
xlim(-0.25,0.25)
show()
|
N=fix(round(T/dtsuper));
fsuper=(mod((0:N-1)+fix(N/2),N)-fix(N/2))/(N*dtsuper);
X=fft(xsuper);
Psuper=abs(X).^2/N^2;
N=fix(round(T/dtabt));
fabt=(mod((0:N-1)+fix(N/2),N)-fix(N/2))/(N*dtabt);
X=fft(xabt);
Pabt=abs(X).^2/N^2;
stem(fsuper,Psuper,'.')
hold on
plot(fabt,Pabt,'o')
hold off
xlim([-0.25 0.25])
|
| korrekte Rekonstruktion des niederfreqenten Anteils |
dtrek=0.1
trek=arange(0,int(round(T/dtrek)))*dtrek
xrek=zeros(len(trek))
for j in range(0,len(xrek)):
if abs(trek[j]/dtabt-round(trek[j]/dtabt))>1e-9:
if len(xabt)%2==0:
for i in range(0,len(xabt)):
xrek[j]+=xabt[i]*sin(pi*(trek[j]-tabt[i])/dtabt)*cos(pi*(trek[j]-tabt[i])/dtabt/float(len(xabt)))/sin(pi*(trek[j]-tabt[i])/dtabt/float(len(xabt)))/float(len(xabt))
else:
for i in range(0,len(xabt)):
xrek[j]+=xabt[i]*sin(pi*(trek[j]-tabt[i])/dtabt)/sin(pi*(trek[j]-tabt[i])/dtabt/float(len(xabt)))/float(len(xabt))
else:
xrek[j]=xabt[int(round(trek[j]/dtabt))%len(xabt)]
plot(tsuper,xsuper,tsuper,xfilt,tabt,xabt,'o')
plot(trek,xrek,'--',dashes=[10,10])
show()
|
dtrek=0.1;
trek=(0:fix(round(T/dtrek))-1)*dtrek;
xrek=zeros(1,length(trek));
for j=1:length(xrek)
if abs(trek(j)/dtabt-fix(trek(j)/dtabt))>1e-9
if mod(length(xabt),2)==0
for i=1:length(xabt)
xrek(j)=xrek(j)+xabt(i)*sin(pi*(trek(j)-tabt(i))/dtabt)*cos(pi*(trek(j)-tabt(i))/dtabt/length(xabt))/sin(pi*(trek(j)-tabt(i))/dtabt/length(xabt))/length(xabt);
end
else
for i=1:length(xabt)
xrek(j)=xrek(j)+xabt(i)*sin(pi*(trek(j)-tabt(i))/dtabt)/sin(pi*(trek(j)-tabt(i))/dtabt/length(xabt))/length(xabt);
end
end
else
xrek(j)=xabt(mod(fix(trek(j)/dtabt),length(xabt))+1);
end
end
plot(tsuper,xsuper,tsuper,xfilt,tabt,xabt,'o',trek,xrek,'--')
|