fe_curve

SDT-base Contents Functions

fe_curve

Purpose

Generic handling of curves and signal processing utilities

Syntax

  out=fe_curve('command',MODEL,'Name',...);

Commands

fe_curve is used to handle curves and do some basic signal processing. The format for curves is described in section 7.9. The iiplot interface may be used to plot curves and a basic call would be iiplot(Curve) to plot curve data structure Curve.

Accepted commands are

bandpass Unit f_min f_max

out=fe_curve('BandPass Unit f_min f_max',signals);
realizes a true bandpass filtering (i.e. using fft() and ifft()) of time signals contained in curves signals. f_min and f_max are given in units Unit, whether Hertz(Hz) or Radian(Rd). With no Unit, f_min and f_max are assumed to be in Hertz.

% apply a true bandpasss filter to a signal
out=fe_curve('TestFrame');% 3 DOF oscillator response to noisy input
fe_curve('Plot',out{2});  % "unfiltered" response
filt_disp=fe_curve('BandPass Hz 70 90',out{2}); % filtering 
fe_curve('Plot',filt_disp); title('filtered displacement');

datatype [,cell]

out=fe_curve('DataType',DesiredType);
returns a data structure describing the data type, useful to fill .xunit and .yunit fields for curves definition. DesiredType could be a string or a number corresponding to the desired type. With no DesiredType, the current list of available types is displayed. One can specify the unit with out=fe_curve('DataType',DesiredType,'unit');.

DataTypeCell returns a cell array rather than data structure to follow the specification for curve data structures. Command option -label"lab" allows directly providing a custom label named lab in the data type.

getCurve,Id,IdNew

GetCurve: recover curve by name curve=fe_curve('getcurve',model,'curve_name');
extracts curve curve_name from model.Stack or the possible curves attached to a load case. If the user does not specify any name, all the curves are returned in a cell array.
GetIdval: recover curve by .ID field, equal to val. curve=fe_curve('GetId val',model);
GetIdNew: generate a new identifier interger unused in any curve of the model. i1=fe_curve('GetIdNew',model);

h1h2 input_channels

FRF=fe_curve('H1H2 input_channels',frames,'window');
computes H1 and H2 FRF estimators along with the coherence from time signals contained in cell array frames using window window. The time vector is given in frames{1}.X while input_channels tells which columns of in frames{1}.Y are inputs. If more than one input channel is specified, true MIMO FRF estimation is done, and Hν is used instead of H2. When multiple frames are given, a mean estimation of FRF is computed.

Note: To ensure the proper assembly of H1 and Hν in MIMO FRF estimation case, a weighing based on maximum time signals amplitude is used. To use your own, use
FRF=fe_curve('H1H2 input_channels',frames,window,weighing);
where weighing is a vector containing weighing factors for each channel. To avoid weighing, use FRF=fe_curve('H1H2 input_channels',frames,window,0); . For an example see sdtweb('start_time2frf','h1h2')

noise

OBSOLETE : use fe_curve TestNoise instead

noise=fe_curve('Noise',Nw_pt,fs,f_max);
computes a Nw_pt points long time signal corresponding to a "white noise", with sample frequency fs and a unitary power spectrum density until f_max. fs/2 is taken as f_max when not specified. The general shape of noise power spectrum density, extending from 0 to fs/2, can be specified instead of f_max.

% computes a 2 seconds long white noise, 1024 Hz of sampling freq.
% with "rounded" shape PSD    
fs=1024; sample_length=2;
Shape=exp(fe_curve('window 1024 hanning'))-1; 
noise_h=fe_curve('noise',fs*sample_length,fs,Shape);
noise_f=fe_curve('fft',noise_h);
figure(1);
subplot(211);fe_curve('plot -gca',noise_h);axis tight;
subplot(212);fe_curve('plot -gca',noise_f);axis tight;

plot

fe_curve('plot',curve); plots the curve curve.
fe_curve('plot',fig_handle,curve); plots curve in the figure with handle fig_handle.
fe_curve('plot',model,'curve_name'); plots the curve of model.Stack named curve_name.
fe_curve('plot',fig_handle,model,curve_name); plots curve named curve_name stacked in .Stack field of model model.

% Plot a fe_curve signal
% computes a 2 seconds long white noise, 1024 Hz of sampling freq.
fs=1024; sample_length=2;
noise=fe_curve('noise',fs*sample_length,fs);
noise.xunit=fe_curve('DataType','Time');
noise.yunit=fe_curve('DataType','Excit. force');
noise.name='Input force';

fe_curve('Plot',noise);

resspectrum [True, Pseudo] [Abs., Rel.] [Disp., Vel., Acc.]

out=fe_curve('ResSpectrum',signal,freq,damp);
computes the response spectrum associated to the time signal given in signal. Time derivatives can be obtained with option -v or -a. Time integration with option +v or +a. Pseudo derivatives with option PseudoA or PseudoV. freq and damp are frequencies (in Hz) and damping ratios vectors of interest for the response spectra. For example

wd=fileparts(which('d_ubeam'));
% read the acceleration time signal
bagnol_ns=fe_curve(['read' fullfile(wd,'bagnol_ns.cyt')]);

% read reference spectrum
bagnol_ns_rspec_pa= fe_curve(['read' fullfile(wd,'bagnol_ns_rspec_pa.cyt')]);
% compute response spectrum with reference spectrum frequencies
% vector and 5% damping
RespSpec=fe_curve('ResSpectrum PseudoA',...
                  bagnol_ns,bagnol_ns_rspec_pa.X/2/pi,.05);

fe_curve('plot',RespSpec); hold on;
plot(RespSpec.X,bagnol_ns_rspec_pa.Y,'r');
legend('fe\_curve','cyberquake');

returny

If curve has a .Interp field, this interpolation is taken in account. If .Interp field is not present or empty, it uses a degree 2 interpolation by default.

To force a specific interpolation (over passing .interp field, one may insert the -linear, -log or -stair option string in the command.

To extract a curve curve_name and return the values Y corresponding to the input X, the syntax is

y = fe_curve('returny',model,curve_name,X);
Given a curve data structure, to return the values Y corresponding to the input X, the syntax is

y = fe_curve('returny',curve,X);

set

This command sets a curve in the model. 3 types of input are allowed:

A data structure, model=fe_curve(model,'set',curve_name,data_structure)
A string to interprete, model=fe_curve(model,'set',curve_name,string)
A name referring to an existing curve (for load case only), model=fe_curve( model, 'set LoadCurve',load_case,chanel,curve_name). This last behavior is obsolete and should be replaced in your code by a more general call to fe_case SetCurve.

When you want to associate a curve to a load for time integration it is preferable to use a formal definition of the time dependence (if not curve can be interpolated or extrapolated).

The following example illustrates the different calls.

% Sample curve assignment to modal loads in a model
model=fe_time('demo bar'); q0=[];

% curve defined by a by-hand data structure:
c1=struct('ID',1,'X',linspace(0,1e-3,100), ...
     'Y',linspace(0,1e-3,100),'data',[],...
     'xunit',[],'yunit',[],'unit',[],'name','curve 1');
model=fe_curve(model,'set','curve 1',c1);
% curve defined by a string to evaluate (generally test fcn):
model=fe_curve(model,'set','step 1','TestStep t1=1e-3');
% curve defined by a reference curve:
c2=fe_curve('test -ID 100 ricker dt=1e-3 A=1');
model=fe_curve(model,'set','ricker 1',c2);
c3=fe_curve('test eval sin(2*pi*1000*t)'); % 1000 Hz sinus
model=fe_curve(model,'set','sin 1',c3);

% define Load with curve definition
LoadCase=struct('DOF',[1.01;2.01],'def',1e6*eye(2),...
            'curve',{{fe_curve('test ricker dt=2e-3 A=1'),...
                      'ricker 1'}});
model = fe_case(model,'DOFLoad','Point load 1',LoadCase);

% modify a curve in the load case
model=fe_case(model,'SetCurve','Point load 1','TestStep t1=1e-3',2);

% the obsolete but supported call was
model=fe_curve(model,'set LoadCurve','Point load 1',2,'TestStep t1=1e-3');

% one would prefer providing a name to the curve, 
% that will be stacked in the model
model=fe_case(model,'SetCurve','Point load 1',...
 'my\_load','TestStep t1=1e-3',2);

Test ...

The test command handles a large array of analytic and tabular curves. In OpenFEM all parameters of each curve must be given in the proper order. In SDT you can specify only the ones that are not the default using their name.

When the abscissa vector (time, frequency, ...) is given as shown in the example, a tabular result is returned.

Without output argument the curve is simply plotted.

% Standard generation of parametered curves
fe_curve('test')  % lists curently implemented curves

t=linspace(0,3,1024); % Define abscissa vector
% OpenFEM format with all parameters (should be avoid):
C1=fe_curve('test ramp 0.6 2.5 2.3',t);
C2=fe_curve('TestRicker 2 2',t);

% SDT format non default parameters given with their name
%  definition is implicit and will be applied to time vector
%  during the time integration: 
C3=fe_curve('Test CosHan f0=5 n0=3 A=3'); 
C4=fe_curve('testEval 3*cos(2*pi*5*t)');

% Now display result on time vector t:
C3=fe_curve(C3,t);C4=fe_curve(C4,t)
figure(1);plot(t,[C1.Y C2.Y C4.Y C3.Y]);
legend(C1.name,C2.name,C4.name,C3.name)

A partial list of accepted test curves follows

Testsin, Testcos, TestTan, TestExp, accept parameters T period and A amplitude. -stoptime Tf will truncate the signal.
TestRamp t0=t0 t1=t1 Yf=Yf has a ramp starting at zero until t0 and going up to Yf at t1. The number of intermediate value can be controlled with the abscissa vector.
To define a gradual load, for non linear static for example, a specific call with a Nstep parameter can be performed : TestRamp NStep=NStep Yf=Yf. For example, to define a 20 gradual steps to 1e-6 :R1=fe_curve('TestRamp NStep=20 Yf=1e-6');
TestRicker dt=dt A=A t0=t0 generates a Ricker function typically used to represent impacts of duration dt and amplitude A, starting from time t0.
TestSweep fmin=fmin fmax=fmax t0=t0 t1=t1 generates a sweep cosine from t0 to t1, with linear frequency sweeping from f0 to f1.
Y=cos(2*pi*(fmin+(fmax−fmin)*t−t0/t1−t0)*(t−t0)) for t0<t<t1, Y=0 elsewhere.
TestStep t1=t1 generates a step which value is one from time 0 to time t1.
TestNoise -window"window" computes a time signal corresponding to a white noise, with the power spectrum density specified as the window parameter. For example TestNoise "Box A=1 min=0 max=200" defines a unitary power spectrum density from 0 Hz to 200 Hz.
TestBox A=A min=min max=max generates a sample box signal from min to max abscissa, with an amplitude A.
TestEval str generates the signal obtained by evaluating the string str function of t. For example R1=fe_curve('Test eval sin(2*pi*1000*t)',linspace(0,0.005,501)); iiplot(R1)

One can use fe_curve('TestList') to obtain a cell array of the test keywords recognized.

testframe

out=fe_curve('TestFrame'); computes the time response of a 3 DOF oscillator to a white noise and fills the cell array out with noise signal in cell 1 and time response in cell 2. See sdtweb fe_curve('TestFrame') to open the function at this example.

timefreq

out=fe_curve('TimeFreq',Input,xf);
computes response of a system with given transfer functions FRF to time input Input. Sampling frequency and length of time signal Input must be coherent with frequency step and length of given transfer FRF.

% Plot time frequency diagrams of signals
fs=1024; sample_length=2;                   % 2 sec. long white noise
noise=fe_curve('noise',fs*sample_length,fs);% 1024 Hz of sampling freq.
[t,f,N]=fe_curve('getXTime',noise);

% FRF with resonant freq. 50 100 200 Hz, unit amplitude, 2% damping 
xf=nor2xf(2*pi*[50 100 200].',.02,[1 ; 1 ; 1],[1 1 1],2*pi*f);

Resp=fe_curve('TimeFreq',noise,xf); % Response to noisy input
fe_curve('Plot',Resp); title('Time response');

Window ...

Use fe_curve window to list implemented windows. The general calling format is win=fe_curve('Window Nb_pts Type Arg'); which computes a Nb_pts points window. The default is a symmetric window (last point at zero), the command option -per clips the last point of a N+1 long symmetric window.

For the exponential window the arguments are three doubles. win = fe_curve('Window 1024 Exponential 10 20 10'); returns an exponential window with 10 zero points, a 20 point flat top, and a decaying exponential over the 1004 remaining points with a last point at exp(-10).

win = fe_curve('Window 1024 Hanning'); returns a 1024 point long hanning window.