demo10.m

来自「这是一个 York大学开发的FDTD代码。很好用的工具」· M 代码 · 共 54 行

% 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;
title('Sum of damped Sine waves')
xlabel('Time')
plot(t,x)
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);
title('Sum of damped Sine waves extended by method of Prony')
xlabel('Time')
plot(t2,real(x2),t,x)

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); % Original series
XD(1)=0.5*XD(1); 
X2D=2*Ts*fft(x2); %Prony Extended series
X2D(1)=0.5*X2D(1);
x3=[x,zeros(1,N2-N)]; % Series extended with zeros
X3D=2*Ts*fft(x3);  
f2=linspace(0,fs*(N2-1)/N2,N2);

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