matlab: fix DFT for periodic signals

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>
pull/14/head
Thorsten Liebig 2015-12-13 16:27:10 +01:00
parent 93c878b28f
commit 1eca3dd23a
4 changed files with 39 additions and 11 deletions

View File

@ -1,32 +1,47 @@
function f_val = DFT_time2freq( t, val, freq ) function f_val = DFT_time2freq( t, val, freq, signal_type )
% f_val = DFT_time2freq( t, val, freq ) % f_val = DFT_time2freq( t, val, freq, signal_type )
% %
% computes the DFT at the given frequencies % computes the DFT at the given frequencies
%
% parameter:
% t : time vector
% val: data vector
% freq: DFT frequency vector
% signal_type: 'pulse' (default), 'periodic'
%
% return values:
% f_val: single-sided spectrum % f_val: single-sided spectrum
% %
% example: % example:
% t=linspace(0,1,100); % t=linspace(0,1,100);
% t_val=0.9*sin(2*pi*3*t); % sine wave; amplitude 0.9; frequency 3 Hz % t_val=0.9*sin(2*pi*3*t); % sine wave; amplitude 0.9; frequency 3 Hz
% f=linspace(1,5,101); % f=linspace(1,5,101);
% f_val=DFT_time2freq( t, t_val, f ); % f_val=DFT_time2freq( t, t_val, f, 'periodic' );
% interp1(f,abs(f_val),3) % interp1(f,abs(f_val),3)
% ans = 0.8910 % ans = 0.8910
% plot( t, t_val ) % plot( t, t_val )
% plot( f, abs(f_val) ) % plot( f, abs(f_val) )
%
% See also FFT_time2freq
if numel(t) ~= numel(val) if numel(t) ~= numel(val)
error 'numel(t) ~= numel(val)' error 'numel(t) ~= numel(val)'
end end
dt = t(2)-t(1); if nargin<4
signal_type = 'pulse';
end
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
if strcmpi(signal_type, 'pulse')
dt = t(2)-t(1);
f_val = f_val * dt; f_val = f_val * dt;
elseif strcmpi(signal_type, 'periodic')
f_val = f_val / length(t);
else
error 'unknown signal type'
end
f_val = f_val * 2; % single-sided spectrum f_val = f_val * 2; % single-sided spectrum

View File

@ -1,6 +1,8 @@
function [f,val] = FFT_time2freq( t, val ) function [f,val] = FFT_time2freq( t, val )
% [f,val] = FFT_time2freq( t, val ) % [f,val] = FFT_time2freq( t, val )
% %
% Note: This function can only be used for pulse signals
%
% See also DFT_time2freq % See also DFT_time2freq
dt=t(2)-t(1); % timestep dt=t(2)-t(1); % timestep

View File

@ -36,15 +36,22 @@ if (nargin<2)
end end
AR_order = 0; AR_order = 0;
SignalType = 'pulse';
for n=1:2:numel(varargin) for n=1:2:numel(varargin)
if (strcmp(varargin{n},'AR')==1) if (strcmp(varargin{n},'AR')==1)
AR_order = varargin{n+1}; AR_order = varargin{n+1};
elseif strcmpi(varargin{n},'SignalType')
SignalType = varargin{n+1};
else else
warning('CSXCAD:ReadUI', ['"' varargin{n} '" is an unknown argument']); warning('CSXCAD:ReadUI', ['"' varargin{n} '" is an unknown argument']);
end end
end end
if strcmpi(SignalType,'periodic') && AR_order>0
error 'auto-regressive model not compatible with periodic signals'
end
if (ischar(files)) if (ischar(files))
filenames{1}=files; filenames{1}=files;
else else
@ -66,6 +73,9 @@ for n=1:numel(filenames)
end end
if (nargin<3) || isempty(freq) if (nargin<3) || isempty(freq)
if strcmpi(SignalType,'periodic')
warning 'ReadUI: periodic signal type not supported by FFT'
end
[UI.FD{n}.f,UI.FD{n}.val] = FFT_time2freq( t,val ); [UI.FD{n}.f,UI.FD{n}.val] = FFT_time2freq( t,val );
else else
UI.FD{n}.f = freq; UI.FD{n}.f = freq;
@ -82,15 +92,15 @@ for n=1:numel(filenames)
end end
if (EC~=0) if (EC~=0)
warning('CSXCAD:ReadUI','AR estimation failed, skipping...') warning('CSXCAD:ReadUI','AR estimation failed, skipping...')
UI.FD{n}.val = DFT_time2freq( t, val, freq ); UI.FD{n}.val = DFT_time2freq( t, val, freq, SignalType );
end end
elseif (AR_order<=0) elseif (AR_order<=0)
UI.FD{n}.val = DFT_time2freq( t, val, freq ); UI.FD{n}.val = DFT_time2freq( t, val, freq, SignalType );
else else
[val_ar t_ar UI.FD{n}.val EC] = AR_estimate( t, val, freq, AR_order); [val_ar t_ar UI.FD{n}.val EC] = AR_estimate( t, val, freq, AR_order);
if (EC~=0) if (EC~=0)
warning('CSXCAD:ReadUI','AR estimation failed, skipping...') warning('CSXCAD:ReadUI','AR estimation failed, skipping...')
UI.FD{n}.val = DFT_time2freq( t, val, freq ); UI.FD{n}.val = DFT_time2freq( t, val, freq, SignalType );
end end
end end
end end

View File

@ -19,6 +19,7 @@ function [port] = calcPort( port, SimDir, f, varargin)
% 'RefPlaneShift': for transmission lines only, See also calcTLPort for % 'RefPlaneShift': for transmission lines only, See also calcTLPort for
% more details % more details
% 'SwitchDirection': 0/1, switch assumed direction of propagation % 'SwitchDirection': 0/1, switch assumed direction of propagation
% 'SignalType': 'pulse' (default) or 'periodic'
% %
% output: % output:
% % output signals/values in time domain (TD): % % output signals/values in time domain (TD):