Time response - comments are welcome

General usage
1

Today I simulated the time response of a sine burst through a Butterworth filter. This is the first time that I do time-response calculations with Scilab. Comments and criticism is welcome.

// roll_off.sce
//
// Simulate a sine-burst signal through a 4th-order Butterworth high-pass
// filter with a corner-frequency at 50 Hz (e.g., a bass-reflex loudspeaker)
// to illustrate phase distortion.

clc; clear; close(winsid());    // use clf(): ??
saveplot = %f;
plotscale = 1.5;    // plotscale = 1.0 provides default figure size

// Extract the filename and path of this script:
[units,typs,nams]=file(); // nams(2) = script file name incl. path
fpathname=strsplit(nams(2), [filesep()]);
scriptpath = get_absolute_file_path(fpathname($));

sr = 20000; // 20 kHz sample rate (not critical, 2 kHz also works)
fc = 10;    // Hz, cosine window, create a 10 Hz = 0.1 second burst
fs = 50;    // Hz, sine burst center frequency
tau = 0:fc/sr:1-1/sr;   // generate sr samples all < 1 (tau, normalized time)
t = tau/fc; // absolute time in seconds, depending on window size
sb = 0.5*(1 - cos(2*%pi*fc*t)).*sin(2*%pi*fs*t);    // sb = sine burst
t2 = [tau 1+tau]/fc;    // extend the time x 2
sb2 = [sb zeros(sb)];   // add trailing zeros, so we can see the response tail

// Create the transfer function in Scilab
w0 = 2*%pi*fs;      // Define corner frequency w0 at fs = 50 Hz
s=poly(0,'s')/w0;   // Define the Laplace variable s
a1 = sqrt(8) * cos(%pi/8);
a2 = 2 + sqrt(2);   // Polynomial coefficients of 4th-order Butterworth filter
a3 = a1;
// Normalized Butterworth 4th-order high-pass transfer function:
num = s^4;
den = s^4 + a1*s^3 + a2*s^2 + a3*s + 1;
// Convert P to a continuous-time linear transfer function
P = syslin('c',num,den);

// gainplot(P);        // Show the frequency response, verification

[y, x] = csim(sb2,t2,P);
timeplot = scf();
timeplot.figure_size = plotscale * timeplot.figure_size;
plot(t2,sb2,'-g');    // show input signal in green
plot([t2',t2'],[y',0*t2'],'-b'); // show output signal in blue
xlabel("t (s)");
ylabel("Amplitude level (relative)");
legend("input","output",1);

if saveplot then 
    fnams = scriptpath + "\" + "burst_response.png";
    xs2png(gcf(),fnams);
end

With kind regards,
Claus

1 Like