zcy
/
openEMS
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

matlab: critical bug fix in dft/fft calculations

pull/1/head
Thorsten Liebig 2011-01-31 10:57:20 +01:00
parent 529cbc1305
commit 15826e910e
3 changed files with 27 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
matlab/DFT_time2freq.m
View File

@ -22,11 +22,13 @@ end
% convert absolute time into relative time % convert absolute time into relative time
t = t - t(1); t = t - t(1);
dt = t(2)-t(1);
f_val = zeros(1,numel(freq)); f_val = zeros(1,numel(freq));
for f_idx=1:numel(freq) for f_idx=1:numel(freq)
f_val(f_idx) = sum( val .* exp( -1i * 2*pi*freq(f_idx) * t ) ); f_val(f_idx) = sum( val .* exp( -1i * 2*pi*freq(f_idx) * t ) );
end end
f_val = f_val / numel(t);
f_val = f_val * dt;
f_val = f_val * 2; % single-sided spectrum f_val = f_val * 2; % single-sided spectrum

2
matlab/FFT_time2freq.m
View File

@ -7,7 +7,7 @@ dt=t(2)-t(1); % timestep
L=numel(val); % signal length L=numel(val); % signal length
NFFT = 2^nextpow2(L); % next power of 2 (makes fft fast) NFFT = 2^nextpow2(L); % next power of 2 (makes fft fast)
%very fine freq resolution... NFFT = NFFT+100000; %very fine freq resolution... NFFT = NFFT+100000;
val = fft( val, NFFT)/L; val = fft( val, NFFT)*dt;
f = 1/(2*dt) * linspace(0,1,NFFT/2+1); f = 1/(2*dt) * linspace(0,1,NFFT/2+1);
val = 2*val(1:NFFT/2+1); % single-sided spectrum val = 2*val(1:NFFT/2+1); % single-sided spectrum

30
matlab/examples/other/gauss_excitation_test.m
View File

@ -4,27 +4,28 @@
clear clear
close all close all
clc
f0 = 0; f0 = 0e9;
fc = 10e9; fc = 10e9;
dT = 8e-12; % sample time-step dT = 1e-12; % sample time-step
sigma = 1/sqrt(8/9)/pi/fc;
t0 = sqrt(18)/sqrt(8/9)/pi/fc;
len = 2 * 9/(2*pi*fc) / dT; % gauss length len = 2 * 9/(2*pi*fc) / dT; % gauss length
for n=1:len for n=1:len
ex(n)=cos(2*pi*f0*((n-1)*dT - 9/(2*pi*fc))) .* exp(-1*(2*pi*fc*(n-1)*dT/3-3).^2); t_(n) = (n-1)*dT;
t_(n)=(n-1)*dT; ex(n) = cos(2*pi*f0*((n-1)*dT - 9/(2*pi*fc))) .* exp(-((t_(n)-t0)/sigma)^2/2);
end end
plot(t_/1e-9,ex) plot(t_/1e-9,ex)
xlabel( 'time (ns)' ); xlabel( 'time (ns)' );
ylabel( 'amplitude' ); ylabel( 'amplitude' );
disp( ['Amplitude at t=0: ' num2str(20*log10(abs(ex(1))/1)) ' dB'] ); disp( ['Amplitude at t=0: ' num2str(20*log10(abs(ex(1))/1)) ' dB'] );
val = DFT_time2freq( t_, ex, [f0-fc f0 f0+fc] ); val = DFT_time2freq( t_, ex, [f0-fc f0 f0+fc] );
@ -43,10 +44,22 @@ val_fft = val_fft((f0-fc<=f) & (f<=f0+fc));
f = f((f0-fc<=f) & (f<=f0+fc)); f = f((f0-fc<=f) & (f<=f0+fc));
hold on hold on
plot( f/1e9, abs(val_fft), 'r' ) plot( f/1e9, abs(val_fft), 'r' )
hold on
if (f0==0)
Fw = sigma*sqrt(2*pi)*exp(-0.5*(sigma*2*pi*f).^2);
plot( f/1e9, 2*abs(Fw), 'g--' )
legend('dft','fft','analytic')
else
legend('dft','fft')
end
xlim([0 max(f)/1e9])
xlabel( 'frequency (GHz)' ); xlabel( 'frequency (GHz)' );
ylabel( 'amplitude' ); ylabel( 'amplitude' );
% dB % dB
figure figure
val = val(freq>=0); val = val(freq>=0);
@ -54,3 +67,6 @@ freq = freq(freq>=0);
plot( freq/1e9, 20*log10(abs(val)/max(abs(val))), 'r' ) plot( freq/1e9, 20*log10(abs(val)/max(abs(val))), 'r' )
xlabel( 'frequency (GHz)' ); xlabel( 'frequency (GHz)' );
ylabel( 'amplitude (dB)' ); ylabel( 'amplitude (dB)' );