[Scilab-users] Matlab to Scilab

[Claus Futtrup]

Hi there

A colleague of mine cannot write Scilab code yet, so he writes in Matlab 
and I convert. So far I've run into a couple of issues that I don't know 
what to make of. Please help us. The .sce script file is attached.

1) The "plot" which is commented out on line 41 crashes Scilab. I wonder 
why - is it because I'm plotting 24000 data points?

2) At the end of the script (line 111) I get a memory allocation error 
with "conv" ... please help on this one.

P.S. My colleague is trying to generate IEC 268-1 noise ... later to 
compress it and make IEC 268-5 (EN 60268-5) test signal for power 
testing loudspeakers.

/Claus


---
This email is free from viruses and malware because avast! Antivirus protection is active.
http://www.avast.com
-------------- next part --------------
// IEC noise generator and filter
// used to calculate voltages during power tests
// creater: PB + CF
xdel(winsid()); // Clear graphics
// clc
clear; // Clear variables
//close all
tic

// define state variables
Fs = 48000; // sampling freq
T = 60; // length of signal in seconds

// Create sample array
t = linspace(1,Fs*T,Fs*T);

// generate white gausian noise
rand_array = grand(1,length(t),"nor",1,1);
// rand_array = zeros(1,length(t)); rand_array(1) = 1; %dirac impulse

// specify frequency weighting 20Hz - 20kHz
// n = -17:1:13;
f_weighting = [0       20    25   31.5   40   50   63   80   100  125  160  200 250 315 400 500 630 800 1000 1250 1600 2000 2500 3150 4000 5000 6300 8000 10000 12500 16000 20000 Fs/2]; // 1000 * (10.^(n/10));
weighting =   [-100   -13.5 -10.2 -7.4  -5.2 -3.5 -2.3 -1.4 -0.9 -0.5 -0.2 -0.1 0   0   0   0   0   0   -0.1 -0.3 -0.6 -1.0 -1.6 -2.5 -3.7 -5.1 -7.0 -9.4 -11.9 -14.8 -18.2 -21.6 -100];
limits =      [%inf      3     2     1     1    1    1    1    0.8  0.6  0.5  0.5 0.5 0.5 0.5 0.5 0.5 0.5  0.6  0.7  0.8  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0   1.5   2.0   3.0  %inf];

a = gca(); // prepare plot of figure a
plot(f_weighting(2:$-1),weighting(2:$-1)); // plot graph a - of weighting
a.log_flags = "lnn"; // make x-axis logarithmic
legend("Weighting curve"); // add curve description to legend
// a.data_bounds = []; // limit plot range
// grid on - there's no Scilab equivalent, AFAIK.

// interpolate weighting array to get more frequency points
N = Fs; // number of points
f_weighting_N = linspace(min(f_weighting),max(f_weighting),N/2);
weighting_N = interp1(f_weighting,weighting,f_weighting_N);
weighting_N(1) = weighting(1);

// b = gca();
// plot(f_weighting_N, weighting_N);
// b.log_flags = "lnn";

// calculate magnitude of weighting in linear scale
IEC_noise_weighting = 10.^ (weighting_N./20);

// create freq array for IFFT
W = zeros(1,N);
inv_IEC_noise_weighting = conj(IEC_noise_weighting(:,$:-1:1));
W(1:N/2) = IEC_noise_weighting;
W(N/2 + 1:N) = inv_IEC_noise_weighting;

// get impulse response of IEC weighting filter
T_weighting = real(ifft(W));

// shift the impulse response back in time to make it causal
IEC_filter_impulse = zeros(1,N);
IEC_filter_impulse(1:(N/2)) = T_weighting((N/2)+1:N);
IEC_filter_impulse((N/2)+1:N) = T_weighting(1:(N/2));

function [B] = rot90(A,k)
    // define rot90 as a function call (not native) in Scilab

  // Output variables initialisation (not found in input variables)
  B=[];

  // Number of arguments in function call
  [%nargout,%nargin] = argn(0)

  // Display mode
  mode(0);

  // Display warning for floating point exception
  ieee(1);

  if ~type(A)==1 then
    error("rot90: Wrong type for input argument #1: Real or complex matrix expected.");
    return;
  end;
  [m,n] = size(A);
  if %nargin==1 then
    k = 1;
  else
    if ~type(k)==1 then
      error("rot90: Wrong type for input argument #2: A real expected.");
      return;
    end;
    k = k-fix(k ./4) .*4;
    if mtlb_logic(k,"<",0) then
      k = mtlb_a(k,4);
    end;
  end;

  if mtlb_logic(k,"==",1) then
    A = mtlb_0(A);
    B = A(n:-1:1,:);
  elseif mtlb_logic(k,"==",2) then
    B = A(m:-1:1,n:-1:1);
  elseif mtlb_logic(k,"==",3) then
    B = A(m:-1:1,:);
    B = mtlb_0(B);
  else
    B = A;
  end;
endfunction

// window the impulse response
IEC_filter_impulse_W = rot90(window('hn',N)) .* IEC_filter_impulse';

// convolve IEC response filter with white noise
IEC_noise = conv(rand_array,IEC_filter_impulse_W);
D = length(IEC_noise); // length of array

// find RMS level of IEC noise
// simple method
RMS_level = sqrt((1/D).*sum(IEC_noise.^2));
// integral method
RMS_int = sqrt((1/D).*trapz(IEC_noise.^2));

// scale IEC noise to unity RMS level
IEC_noise_scaled = IEC_noise ./ RMS_level; // block for impulse signal
// IEC_noise_scaled = IEC_noise; %block for noise signal - use for impulse