demo10.m

这是一个 York大学开发的FDTD代码。很好用的工具

% uses setup from demo0
disp('Prony extension of decaying cosine-waves')

disp('Generate time and frequency steps')
Ts=Tp/(N)
fs=1/Ts
f=linspace(0,fs*(N-1)/N,N);
disp('Generate time steps for one period (t)');
t=linspace(0,Tp*(N-1)/N,N);



disp('Generate time sequence (x) and plot it against time (t)')
disp('Frequencies:')
f0
disp('Time constants:')
Tor
x=cos(2*pi*f0'*t).*exp(-(1./Tor')*t);
x=ones(1,size(f0,2))*x;
plot(t,x)
title('Sum of damped Sine waves')
xlabel('Time')
input('Any key to continue')

disp('Use method of prony to generate/extend time series')
Psiz=input('Number of Prony coeffs? ')
[C,Sp]=prony(x',Ts,Psiz);

disp('Current length='),N
N2=input('length of extended sequence ? ')
x2=invprony(C,Sp,Ts,1,N2);
t2=linspace(0,Tp*(N2-1)/N,N2);

plot(t,x,t2,real(x2))
title('Sum of damped Sine waves extended by method of Prony')
xlabel('Time')
input('Any key to continue')

disp('Generate DFT (X) of time series (x) and Compare with Fourier transform of continous data');
XD=2*Ts*fft(x);
XD(1)=0.5*XD(1);
X2D=2*Ts*fft(x2);
X2D(1)=0.5*X2D(1);
x3=[x,zeros(1,N2-N)]; %Generate truncated FFT too !
X3D=2*Ts*fft(x3);  
f2=linspace(0,fs*(N2-1)/N2,N2);

plot(f(1:N/2)/fs,db(XD(1:N/2)),'o',f0/fs,db(ones(1,size(f0,2)).*Tor),'+',f2(1:N2/2)/fs,db(X2D(1:N2/2)),'-',f2(1:N2/2)/fs,db(X3D(1:N2/2)),'.')
title('Discrete Fourier Transform Spectra (truncation effects)')
xlabel('Frequency/fs')
ylabel('abs(X)')
disp('FFT of: o-Original time series, +-Cosine freqs, --Prony,  .-Zero padded')