[Scilab-users] Frequency response

[Claus Futtrup]

Hi Rafael

My problem is defined as a response function, where I use contour 
integral to achieve the step response. Then I differentiate it to 
achieve the impulse response. In essense I already "know" my response 
function. It contains some special features (it could for example be a 
logarithmic term) and e.g. "syslin" or "bode" will choke on it. If I 
already know the response function, why bother with the FFT? Well, I'd 
like to ensure the time response and frequency response are correct 
(without errors on my part) and one way of verifying is to 
close-the-loop and calculate the FFT, see that input = output. I'd like 
to verify that the absolute levels are correct.

Attached please find the script (10 kb). It's fully functional. 
Apologies for the comments in the script, I didn't clean it up for 
publishing...

Best regards,
Claus

On 15.09.2018 22:47, Rafael Guerra wrote:
> Forgetting MATLAB for a second, can you define your problem ?
>
> How is your specific frequency response defined? In Laplace, Fourier, 
> z-transform, ..., domain?
> ------------------------------------------------------------------------
> *From:* users <users-bounces at lists.scilab.org> on behalf of Claus 
> Futtrup <cfuttrup at gmail.com>
> *Sent:* Saturday, September 15, 2018 11:15:08 PM
> *To:* users at lists.scilab.org
> *Subject:* Re: [Scilab-users] Frequency response
> Hi Samuel and Scilabers
>
> My problem with freq and repfreq is that they require a "sys" input, 
> which implies I need to describe a response function of some sort. For 
> example (from the Scilab web-help):
>
> https://help.scilab.org/docs/6.0.1/en_US/freq.html (example code)
> s=poly <https://help.scilab.org/docs/6.0.1/en_US/poly.html>(0,'s');
> sys=(s+1)/(s^3-5*s+4)
> rep=freq(sys("num"),sys("den"),[0,0.9,1.1,2,3,10,20])
> https://help.scilab.org/docs/6.0.1/en_US/repfreq.html (syntax)
> [ [frq,] repf]=repfreq(sys,fmin,fmax[,step])
> [ [frq,] repf]=repfreq(sys[,frq])
> [ frq,repf,splitf]=repfreq(sys,fmin,fmax[,step])
> [ frq,repf,splitf]=repfreq(sys[,frq])
>
> ... So, the thing with the matlab freqz (as the example repeated below 
> shows) is just a basic FFT with sampling, etc:
>
> h = rand(1,64); // impulse response (Matlab source code)
> fs = 1000;
> Nfft = 128;
> [H,F] = freqz(h,1,Nfft,fs);
> semilogx(F,mag2db(abs(H))); grid on;
> xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)');
>
> It may very well be that I just don't understand the help page of 
> repfreq. For example it says about sys: "A siso or simo linear 
> dynamical system, in state space, transfer function or zpk 
> representations, in continuous or discrete time." - al right then, it 
> seems pretty capable, but so, what's "zpk" for example? ... my 
> apologies for finding the Scilab help to be cryptic and the examples 
> insufficient for me to solve my problem, hence I ask if someone can 
> take above matlab code and make it work in Scilab.
>
> Best regards,
> Claus
>
> On 15.09.2018 00:32, Samuel Gougeon wrote:
>> Le 14/09/2018 à 20:57, Claus Futtrup a écrit :
>>>
>>> Dear Scilabers
>>>
>>> I have calculated an impulse response and wish to do an FFT to 
>>> achieve the frequency response. I know what to expect. In the matlab 
>>> forum someone asked the same question and was recommended to use 
>>> freqz ... I wonder what would be the equivalent function in Scilab?
>>>
>>> https://www.mathworks.com/matlabcentral/answers/350350-how-to-plot-loudspeaker-frequency-response-from-its-impulse-response 
>>>
>>>
>>> For example, to replicate the code snippet (second answer in above 
>>> link), how to do this in Scilab?
>>>
>>> h = rand(1,64); // impulse response (Matlab source code)
>>> fs = 1000;
>>> Nfft = 128;
>>> [H,F] = freqz(h,1,Nfft,fs);
>>
>> Did you have a look around freq() or repfreq()?
>>
>> We have somewhat the equivalence invfreqz(H,F,m,n,W) <=> 
>> frfit(F*2*%pi, H, n, W) // Scilab
>>
>> So you may look for the reciprocal of Scilab's frfit()
>>
>> HTH
>> Samuel
>>
>>
>>
>> _______________________________________________
>> users mailing list
>> users at lists.scilab.org
>> http://lists.scilab.org/mailman/listinfo/users
>
>
>
> <https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient> 
> 	Virus-free. www.avast.com 
> <https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient> 
>
>
>
>
> _______________________________________________
> users mailing list
> users at lists.scilab.org
> http://lists.scilab.org/mailman/listinfo/users




---
This email has been checked for viruses by Avast antivirus software.
https://www.avast.com/antivirus
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20180916/9138d9ac/attachment.htm>
-------------- next part --------------
// stepresponse.sce
//
// This script calculates the step response of a specified bass reflex system
// design based on the Weideman Contour method as described in the AES article:
// "A Contour Integral Method for Transient Response Calculation"
// J. Audio Eng. Soc., vol. 66, no. 5, pp. 360-368, (May 2018).
//
// The transducer data are from the following AES article (not the datasheet):
// "An Added-Mass Measurement Technique for Transducer Parameter Estimation"
// (because it includes visco-elastic parameters, classic Knudsen LOG model)
// and therefore this script is using the Knudsen LOG model for simulation.
//
// Script by Claus Futtrup and Jeff Candy
// Date created: 2018-08-05
// Last modified: 2018-09-14

clear; xdel(winsid()); clc; // Clear variables, graphics and console

// ====================
// ---- USER INPUT ----
// ====================

name = "SEAS L16RNX 8 ohm in 10 liter box"; // Included in PLOT title

// --- DEFINE TRANSDUCER PARAMETERS
Re = 5.8; // ohm - DC resistance - world's simplest "inductance model"
Bl = 7.034; // Tm (or N/A)
Mms = 15.13; // gram
Sd = 104; // cm2
// Instead of the classical Thiele-Small parameters of Cms, Rms,
// we need three parameters, named R0, L0 and beta to describe the
// viscoelastic damping and compliance of the mechanical system.
R0 = 32.20; // ohm (damping) - for gbeta = 0, insert Res (T/S parameter)
L0 = 60.24; // mH (compliance) - for gbeta = 0, insert Lces (T/S parameter)
gbeta = 0.059; // Viscoelasticity, if gbeta = 0, then without creep (i.e. T/S)
// Above three parameters describe damping (R0) and compliance (L0, gbeta).
// We recommend that they are found using the dual-added-mass technique.
// Notice the zero ('0') in the names of the parameters, which show that these
// values are valid following Knudsen LOG model and _not_ classical T/S.
// The following model equations can degenerate to the Thiele-Small model just
// by setting gbeta = 0, R0 = Res = Bl^2 / Rms and L0 = Lces = Cms * Bl^2;

// --- DEFINE BOX PARAMETERS
Vb = 10; // liter
Fp = 40; // Hz - port tuning frequency, Fp = 0 for closed box simulation
// Add box damping (?) - most important are the absorption losses.

// ===========================
// ---- END OF USER INPUT ----
// ===========================

// --- CONVERT TO SI-UNITS
Mms = Mms / 1e3; // convert to kg
Sd = Sd / 1e4; // convert to m2
L0 = L0 / 1e3; // convert to Henry
Vb = Vb / 1e3; // convert to m3

// --- SOME CONSTANTS
rho = 1.201; // Density of air, kg/m3
c = 343.7; // Speed of sound in air, m/s
p0 = 2e-5; // 20 micro-Pascal, Reference pressure level (SPL = 0 dB)
r = 1; // Distance 1 meter
Vref = sqrt(8); // 2.83 Volt (step response signal input level)
                // Should this be 4 Volt (peak amplitude) instead ??

// --- CALCULATE FURTHER PARAMETERS
Cm0 = L0 / Bl^2; // Compliance (mechanical)
ws = 1/sqrt(Cm0*Mms);
F0 = ws / (2 * %pi);
  // This is _not_ the resonance frequency of the complete speaker, but
  // only the frequency-independent part (= excluding viscoelasticity).
  // We do it this way because the viscoelastic compliance and damping
  // is added as a separate expression to the response function.
Rm0 = Bl^2 / R0; // Damping (mechanical), frequency-independent part
Qm0 = (Mms/Rm0) * ws;
Rme = Bl^2 / Re; // Mechanical equivalent of electrical (DC) resistance
Qe0 = (Mms/Rme) * ws;
Qt0 = (1/Qm0 + 1/Qe0)^-1;
Vas = Cm0 * Sd^2 * rho * c^2; // Volume of Air Equivalent to Cm0
lambda = gbeta * log(10) / (1-gbeta * log(ws)); // Knudsen LOG lambda
                                                // Not used, just curious

// --- NORMALIZE DATA
alpha = Vas / Vb;
h = Fp / F0;
// NOTICE: When using the LOG model, these normalizations
// are in physical size quite different from Thiele/Small.

// --- RESPONSE FUNCTION

function y=resp(s)
    cs = 1 - gbeta * log( s ); // Knudsen LOG viscoelastic compliance
    y = s.^4 ./ ( (s.^2 + h^2) .* (1.0./cs + s/Qt0 + s.^2) + alpha * s.^2);
endfunction
// NOTICE: This defines a 4th order system (bass reflex), utilizing the Knudsen
// LOG model with a viscoelastic compliance (i.e. the output is complex and
// gives both stiffness and damping to the final response). The parameters, h,
// Qt0, alpha and gbeta (greek-beta) are defined with a normalized frequency
// given by the LOG model (i.e. not the resonance frequency of the complete
// speaker). The only difference from the original Knudsen LOG article is that
// we use the natural logarithm (not base-10 logarithm).

// --- WEIDEMAN CONTOUR INTEGRAL

N0 = 12; // resolution
steps = 8; // steps per unit time,
t_max = 16; // time (dimensionless) - t_max = 16 gives 128 data points for FFT
t_step = 1/steps;                  // because steps * t_max = 8 * 16 = 128
t = t_step:t_step:t_max;           // which means frequency-step = 14.6 Hz

mu_c = max([1 h]);
t_c = %pi*N0/(12*mu_c);
nt = length(t);

for i=1:nt do
    ti = t(i);       // Selection defined in
    if ti < t_c then // Eqs (31) and (32)
        mu = %pi*N0/(12*ti);
        N = N0;
        d = 3/N;
    else
        mu = mu_c;
        N = ceil(N0*ti/t_c);
        d = 3/N;
    end

    xsk = 0.0; // Perform sum in Eq (24)
    for k=-N:N do
        u = k*d;
        s = mu * (%i*u + 1)^2;
        dsdu = 2.0*mu*(%i-u);
        xsk = xsk + imag(exp(s*ti) * resp(s)/s * dsdu);
    end
    xsk = xsk * d / (2*%pi);
    xs(i) = xsk; // Collect x-values (the normalized pressure response)
end

// --- UNDO NORMALIZATION
tau = t / ws; // in seconds (x-axis)
p1 = rho * Sd / (2 * %pi * r) * Vref * Bl / (Mms * Re);
p = p1*xs; // in Pascal (y-axis)
tau = [0 tau]; // add startpoint (row vector)
p = [p1; p];   // add startpoint (column vector)

// --- PLOT STEP RESPONSE
g_tr = scf();
plot(tau',p,"-xr"); // Transpose X to prevent WARNING
xgrid(color("grey70"));
title_string = "Step (air pressure) time response, " + name;
title(title_string);
xlabel("Time (seconds)");
ylabel("Pressure (Pascal)");

// should I make a real GUI for this, with a table on the left and the graphs
// on the right ??
// http://www.openeering.com/sites/default/files/LHY_Tutorial_Gui.pdf

// Should I make an FFT to see the SPL-response?
// Sample_frequency = 4096 Hz (for example = 100 * Fs)
// Lower frequency limit 10 Hz - and use 10 Hz bins / steps - or too coarse?
// Create sampled input vector, convert step response to impulse response
// Run y = FFT(V)
// Plot(f,y)
// Or does it make more sense to utilize the response function directly?
// - it seems BODE cannot plot functions with LOG() in it ...

// --- CALCULATE AND PLOT IMPULSE RESPONSE
// - take the derivative to obtain the impulse response
tau_step = t_step / ws; // convert normalized time steps to time (seconds)
  // tau_step is our discrete differentiation step (and freq. resp resolution)
H = (p(2:$) - p(1:$-1)) / tau_step; // impulse response
t = tau(2:$) - tau_step/2; // redefine time intervals for impulse response
    // t is not used. It's for plotting the impulse response, if desired
scf();
plot(t',H); // Transpose X to prevent WARNING
xgrid(color("grey70"));
title_string = "Impulse (air pressure) time response, " + name;
title(title_string);
xlabel("Time (seconds)");
ylabel("Pressure (Pascal)");

// The time scale (t) is now non-zero and a zero cannot be added simply by
// extrapolation. Should I instead interpolate and re-instate the old tau
// time steps, so that extrapolation to zero can be done?
// Is it this time thing that f*cks up my FFT?

// --- FREQUENCY RESPONSE
f = 0.1:0.01:2; // This is needed only for frequency response plot, normalized
w = 2 * %pi * f; // Omega, normalized
iw = %i * w; // like s, but 's' variable name is already used

Y = resp(iw);
MAG = 20 * log10( abs(Y) );
PHS = 180 / %pi * atan(imag(Y) , real(Y));

for i=1:length(Y) // unwrap phase
   if PHS(i) < 0.0 then
      PHS(i) = PHS(i)+360.0 ;
   end
end

n0 = rho * (Sd * Bl)^2 / (Mms^2 * 2 * %pi * c * Re);
SPL_ref = 10 * log10(rho * c / (2 *%pi * p0^2));
SPL = SPL_ref + 10 * log10(n0); // SPL = 20 * log10(p1 / p0);
// mprintf("SPL = %3.2e dB\n",SPL); // Result = 85.7 dB

scf();
subplot(2,1,1);
plot(F0*w,MAG + SPL,'-b');
fm = gca(); // Get current axis
fm.log_flags = "lnn"; // Set x-axis = log
xgrid(color("grey70"));
title_string = "Frequency response, " + name;
title(title_string);
xlabel("Frequency (Hz)");
ylabel("Magnitude (dB)");

subplot(2,1,2);
plot(F0*w,PHS,'-b');
fp = gca(); // Get current axis
fp.log_flags = "lnn"; // Set x-axis = log
xgrid(color("grey70"));
// title_string = "Frequency response, " + name;
// title(title_string);
xlabel("Frequency (Hz)");
ylabel("Phase (deg)");

// - do an FFT of the impulse response
// Resp = fft(H);
// Resp_magnitude = 20 * log10(abs(Resp));
// [db,phi] = dbphi(Y); // convert frequency response into magnitude and phase
// f_max = 1/tau_step; // Hz
// f = f_max * (0:length(Y)-1)/length(Y); //associated frequency vector
// scf();
// plot(f'(2:nt/2),abs(db(2:nt/2))); // Transpose X to prevent WARNING
  // Only plot first half of the data (second half = mirror image)
  // Do not plot first value, where the frequency f = 0 Hz.
// fr = gca(); // Get current axis
// fr.log_flags = "lnn"; // Set x-axis = log
// title_string = "Frequency response, " + name;
// title(title_string);
// xlabel("Frequency (Hz)");
// ylabel("Magnitude (dB)");

// Matlab Answers - https://www.mathworks.com/matlabcentral/answers/350350-how-to-plot-loudspeaker-frequency-response-from-its-impulse-response
// h = rand(1,64); % impulse response
// fs = 1000;
// Nfft = 128;
// [H,F] = freqz(h,1,Nfft,fs);
// semilogx(F,mag2db(abs(H))); grid on;
// xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)');

// What is the Scilab equivalent of Matlab 'freqz' ?

// BODE PLOT - NOTICE this doesn't work because 's' cannot cooperate with log():
// s = poly(0, 's');
// num = s^4;
// cs = 1 - gbeta * log( s ); // Knudsen LOG viscoelastic compliance
// den = ( (s^2 + h^2) * (1/cs + s/Qt0 + s^2) + alpha * s^2);
// h = syslin('c', num/den);
// scf(); bode(h, 0.01, 100);