Re: [Scilab-users] Frequency response

[Tan Chin Luh]

Hi Claus,

If I understand you correctly (as mentioned by Rafael "I do not see how 
you can “close the loop”, other that testing your theoretical model 
against experimental data?"


1. you have a system which you have theoretical result (reference model)
2. you have another set of data which you want to compare with the 
reference model


As you're the subject expert in your own code, I will leave it to you as 
it will take time to look into them. :)


Instead, I illustrate something similar (I think) which you might be 
able to adopt into your own code.


Codes below perform following steps:
1. Create a linear model in s domain.
2. Convert the model to z domain.
3. Get the impulse response from the z model.
4. From the impulse response data, find the frequency response (in your 
case, it should be the data H that you obtained somewhere?)

5. Compare the response in 4 with the model response.

rgds,
CL

// Linear System in Time Domain
s=%s;
G  =  syslin('c',  10*s^2  ,  3*s^2  +  s  +  2);

// Discreate Model
Ts  =  0.05;   // Sampling time for discrete model
SS_z  =  cls2dls(tf2ss(G),Ts);
tz  =  [0:Ts  :50]';
u=ones(tz);
y2_step=dsimul(SS_z,u');   //Step response
u=zeros(tz);u(1)=1;
y2_imp=dsimul(SS_z,u');//Impulse response

// Compute Bode plot from Impulse Response
tau_step  =  Ts;
H  =  y2_imp;
fs  =  1/tau_step;

// Following part is the same code as earlier
b  =  poly(H($:-1:1),"z","coeff");
a  =  %z^(length(H)-1);
Gz  =  syslin('d',b,a);
Gz.dt  =  1/fs;
[F,M]=repfreq(Gz,0.001,fs/2,0.01);

// Comparison
bode(G);  // Ideal frequency response
plot(F(2:10:$),20*log10(M(2:10:$)),'r.');  // generated frequency response with 
repfreq




___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

Re: [Scilab-users] Frequency response

[Claus Futtrup]

Hi Tan, et al.

Thank you for the complete response to my initial request.

This is very helpful for me to understand ... but if this is indeed the 
equivalent of the matlab script, I wonder about if the help I found in 
the matlab forum is helpful or not. The reason I'm still wondering is 
because the matlab forum question is exactly like mine - I wish to find 
the fft of a loudspeaker impulse response and I expect a loudspeaker 
frequency response with roll-off at low frequencies, but now that I see 
the result of the below script with random data for the impulse response 
... I see that the generated output is a "filter" - which is not at all 
what I expected and I think the similar response in the matlab forum was 
also not helpful (I do not have access to matlab / cannot test and or 
verify). Sorry for this. The entire concept may not be useful in my 
case, but Scilab supports advanced FFT routines (FFTW) and it should be 
possible to do what I need.


Rafael is right that the previously attached stepresponse script 
explains in detail my situation. If I take the "H" vector (the impulse 
response) in the script and apply a basic fft(H), it looks like this:


Added script code:

Resp  =  fft(H);
[db,phi]  =  dbphi(Resp);  // convert frequency response into magnitude and 
phase
f_max  =  1/tau_step;  // Hz
f  =  f_max  *  (0:length(H)-1)/length(H);  //associated frequency vector
scf();
plot(f'(2:nt/2),db(2:nt/2));  // Transpose X to prevent WARNING
re  =  gca();
re.log_flags  =  "lnn";
xgrid(color("grey70"));
xlabel("Frequency (Hz)");
ylabel("Magnitude (dB)");

RESULT :

This is not exactly what I had in mind, the "shape" of the curve isn't 
right, but then again - there are some inner-workings I don't know about 
(like, is the dbphi function doing what I think it does?). Maybe I 
cannot use dbphi(), or maybe I just have an error somewhere in the 
calculation of the impulse response. I expected a level of approx. 86 dB 
with a roll-off towards the lower frequencies. Here's the expected 
response (from simulation based on the input response function, i.e. not 
based on impulse response and not by fft):




Best regards,
Claus

On 18.09.2018 04:25, Tan Chin Luh wrote:

On 18/9/2018 5:15 AM, Rafael Guerra wrote:

Chin Luh:  good to know but how does that solve Claus Futtrup specific problem?
___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

just notice that the answer similar as the one Tim provided, just a 
bit confused on Claus' comment on the "random" data:


On 17/9/2018 1:13 AM, Claus Futtrup wrote:
What I like about the Matlab example is that random data is generated 
to represent the impulse response, so this represents "any data" ... 
I need that. If Scilab cannot do it, it's OK. 


Do you refer this "random" data to the example provided?

h = rand(1,64); % impulse response <-- This?
fs  =  1000;
Nfft  =  128;
[H,F]  =  freqz(h,1,Nfft,fs);
plot(F,H);  grid  on;

If so, the equivalent scilab code would be as below:

h  =  rand(1,64);
fs  =  1000;
b  =  poly(h($:-1:1),"z","coeff");
a  =  %z^(length(h)-1);
Gz  =  syslin('d',b,a);
Gz.dt  =  1/fs;
[F,H]=repfreq(Gz,0,500,0.01);
plot(F,H); Do note that I replace the last part of semilogx plot and 
replace with plot for more simple codes.

--
Tan Chin Luh
Trity Technologies Sdn Bhd
Tel : +603 80637737
HP : +6013 3691728


___
users mailing list
users@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
___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

Re: [Scilab-users] Frequency response

[Rafael Guerra]

“just notice that the answer similar as the one Tim provided, just a bit 
confused on Claus' comment on the "random" data”

Claus problem is defined in the script: “stepresponse - 2018-09-16.sce” that 
was forwarded earlier


___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

Re: [Scilab-users] Frequency response

[Tan Chin Luh]

On 18/9/2018 5:15 AM, Rafael Guerra wrote:

Chin Luh:  good to know but how does that solve Claus Futtrup specific problem?
___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

just notice that the answer similar as the one Tim provided, just a bit 
confused on Claus' comment on the "random" data:


On 17/9/2018 1:13 AM, Claus Futtrup wrote:
What I like about the Matlab example is that random data is generated 
to represent the impulse response, so this represents "any data" ... I 
need that. If Scilab cannot do it, it's OK. 


Do you refer this "random" data to the example provided?

h = rand(1,64); % impulse response <-- This?
fs  =  1000;
Nfft  =  128;
[H,F]  =  freqz(h,1,Nfft,fs);
plot(F,H);  grid  on;

If so, the equivalent scilab code would be as below:

h  =  rand(1,64);
fs  =  1000;
b  =  poly(h($:-1:1),"z","coeff");
a  =  %z^(length(h)-1);
Gz  =  syslin('d',b,a);
Gz.dt  =  1/fs;
[F,H]=repfreq(Gz,0,500,0.01);
plot(F,H); Do note that I replace the last part of semilogx plot and 
replace with plot for more simple codes.


--
Tan Chin Luh
Trity Technologies Sdn Bhd
Tel : +603 80637737
HP : +6013 3691728

___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

Re: [Scilab-users] Frequency response

[Rafael Guerra]

Chin Luh:  good to know but how does that solve Claus Futtrup specific problem?
___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

Re: [Scilab-users] Frequency response

[Tan Chin Luh]

Hi,

I was facing similar issue sometime ago while trying to translate from 
matlab to scilab, as both of them have diff way of representation.


In Matlab/Octave, when we are using polynomial, we are defining them in 
vector/matrix form such as : num = [1 2 3 4]. The catch is here: when we 
use the s-domain, it represent 1*s^3 + 2*s^2+3*s+4, while if it is used 
in the function in z-domain, it will become 1 + 2*z^-1 + 3*z^-2 + 4*z^-3.


while in Scilab, the polynomial is represented in polynomial datatype.

So, consider following simple example:
h = [1 2 3 4];   % this is coefficient in 1 + 2*z^-1 + 3*z^-2 + 4*z^-3
fs = 1000;
Nfft = 128;
[H,F] = freqz(h,1,Nfft,fs);
plot(F,H);

In Scilab

h = [1 2 3 4];
fs = 1000;
b = poly(h($:-1:1),"z","coeff");
a = %z^(length(h)-1);
Gz = syslin('d',b,a);
Gz.dt = 1/fs;
[F,H]=repfreq(Gz,0,500,0.01);
plot(F,H);

in line 3 - b = poly(h($:-1:1),"z","coeff");  we create b in poly using 
the coefficients, we have to take note that:

--> order of h must be reversed, as poly take the coeff in ascending order
--> this line will create b as 4 +3*z +2*z^2  +z^3

To match the z-domain format, our a would not be 1 anymore, but z^3 as 
in line 4 of the code above

--> (4 +3*z +2*z^2  +z^3) / z^3  --> 1 + 2*z^-1 + 3*z^-2 + 4*z^-3

After this, the remaining codes should be straight forward.

Hope this helps.

p/s: help zpk will help u on the definition i guess.

Thanks.

rgds,
Chin Luh
___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

Re: [Scilab-users] Frequency response

[Rafael Guerra]

Hi Claus,

In your script, the response function resp(s) is definitely not a simple LTI 
system, not even a ratio of polynomials, because of the logarithmic term.
I am not an expert but I suspect that such highly non-linear transfer functions 
are not handled by the standard Scilab DSP tools such as syslin, etc.

I do not see how you can "close the loop", other that testing your theoretical 
model against experimental data?

PS:
ZPK Zero pole gain system representation: zeros Z, poles P, and gain(s) K.

Regards,
Rafael

From: users  On Behalf Of Claus Futtrup
Sent: Sunday, September 16, 2018 5:49 PM
To: users@lists.scilab.org
Subject: Re: [Scilab-users] Frequency response

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 
<mailto:users-boun...@lists.scilab.org> on 
behalf of Claus Futtrup <mailto:cfutt...@gmail.com>
Sent: Saturday, September 15, 2018 11:15:08 PM
To: users@lists.scilab.org<mailto:users@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@lists.scilab.org<mailto:users@lists.scilab.org>

http://lists.scilab.org/mailman/listinfo/users



[https://ipmcdn.avast.com/images/icons/icon-envelope-tick-round-orange-animated-no-repeat-v1.gif]<https://www.avast.com/sig-email?utm_medium=email_source=link_campaign=sig-email_content=emailclient

Re: [Scilab-users] Frequency response

[Claus Futtrup]

Hi Tim

ZPK can mean many things (Google Search, I suppose) ... but after search 
I understand it's something with Zero Pole ... don't know what the 
letter K stands for, though.


What I like about the Matlab example is that random data is generated to 
represent the impulse response, so this represents "any data" ... I need 
that. If Scilab cannot do it, it's OK.


/Claus

On 16.09.2018 19:01, Tim Wescott wrote:

I didn't answer about ZPK because I didn't know either!

It's not so much a Scilab thing as -- are you getting the signal
processing right?

On Sun, 2018-09-16 at 18:15 +0200, Claus Futtrup wrote:

Hi Tim

  >So, this is complicated.

I admitted from the very beginning, it's probably just me that
doesn't
know how to read the Scilab manual. It cannot be complicated for
someone
who knows and/or understand the manual.

BTW, I also notice that nobody answered the question what ZPK means?
...
ZPK = Zagrepa Plivaci Klub ... at least the Scilab could spell out
what
it stands for.

/Claus

On 15.09.2018 22:50, Tim Wescott wrote:

So, this is complicated.

First, if you have an impulse response, in the form of a vector of
samples, then you can turn it into a FIR filter, and you can find
the
frequency response of that with any of the available Scilab tools.
  I'm
not sure if there's a more direct way, but if you had an impulse
response h = [1, 1/2, 1/4, 1/8], then you can turn that into H = 1
+
0.5*%z^(-1) + 0.25*%z^(-2) + 0.125*%z^(-3) (note that %z is a
built-in:
%z = poly(0, 'z')).  Then you can use techniques surrounding
repfreq to
get the frequency response.

Assign a sampling rate to it (i.e., by setting H.dt = 0.001, for a
1kHz
sampling rate) and you can just use the Bode plot function:
bode(H).

Second, you can take your impulse response, pad it with zeros, and
then
take the FFT of it -- i.e. H = fft([h zeros(1, 252)]).  You need to
pad
it because the FFT is only exact for a periodic function in time --
padding it gives you an approximation of an actual impulse response
(and before someone jumps in and says you need to window it -- no,
you
don't, windowing is for sample vectors of continuous signals, not
impulse responses).  But if you don't know WHY the method is
approximate, you're going to have trouble correctly massaging the
data
before it's presented to the FFT, and with interpreting the
results.

On Sat, 2018-09-15 at 22:15 +0200, Claus Futtrup wrote:

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(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
-plo
t-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@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

Virus-free. www.avast.com
   ___

Re: [Scilab-users] Frequency response

[Tim Wescott]

I didn't answer about ZPK because I didn't know either!

It's not so much a Scilab thing as -- are you getting the signal
processing right?

On Sun, 2018-09-16 at 18:15 +0200, Claus Futtrup wrote:
> Hi Tim
> 
>  >So, this is complicated.
> 
> I admitted from the very beginning, it's probably just me that
> doesn't 
> know how to read the Scilab manual. It cannot be complicated for
> someone 
> who knows and/or understand the manual.
> 
> BTW, I also notice that nobody answered the question what ZPK means?
> ... 
> ZPK = Zagrepa Plivaci Klub ... at least the Scilab could spell out
> what 
> it stands for.
> 
> /Claus
> 
> On 15.09.2018 22:50, Tim Wescott wrote:
> > 
> > So, this is complicated.
> > 
> > First, if you have an impulse response, in the form of a vector of
> > samples, then you can turn it into a FIR filter, and you can find
> > the
> > frequency response of that with any of the available Scilab tools.
> >  I'm
> > not sure if there's a more direct way, but if you had an impulse
> > response h = [1, 1/2, 1/4, 1/8], then you can turn that into H = 1
> > +
> > 0.5*%z^(-1) + 0.25*%z^(-2) + 0.125*%z^(-3) (note that %z is a
> > built-in:
> > %z = poly(0, 'z')).  Then you can use techniques surrounding
> > repfreq to
> > get the frequency response.
> > 
> > Assign a sampling rate to it (i.e., by setting H.dt = 0.001, for a
> > 1kHz
> > sampling rate) and you can just use the Bode plot function:
> > bode(H).
> > 
> > Second, you can take your impulse response, pad it with zeros, and
> > then
> > take the FFT of it -- i.e. H = fft([h zeros(1, 252)]).  You need to
> > pad
> > it because the FFT is only exact for a periodic function in time --
> > padding it gives you an approximation of an actual impulse response
> > (and before someone jumps in and says you need to window it -- no,
> > you
> > don't, windowing is for sample vectors of continuous signals, not
> > impulse responses).  But if you don't know WHY the method is
> > approximate, you're going to have trouble correctly massaging the
> > data
> > before it's presented to the FFT, and with interpreting the
> > results.
> > 
> > On Sat, 2018-09-15 at 22:15 +0200, Claus Futtrup wrote:
> > > 
> > > 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(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
> > > > > -plo
> > > > > t-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
> 

Re: [Scilab-users] Frequency response

[Alan Teeder]

-Original Message- 
From: Claus Futtrup 
Sent: Sunday, September 16, 2018 5:15 PM 
To: users@lists.scilab.org 
Subject: Re: [Scilab-users] Frequency response 


Hi Tim


So, this is complicated.


I admitted from the very beginning, it's probably just me that doesn't 
know how to read the Scilab manual. It cannot be complicated for someone 
who knows and/or understand the manual.


BTW, I also notice that nobody answered the question what ZPK means? ... 
ZPK = Zagrepa Plivaci Klub ... at least the Scilab could spell out what 
it stands for.


/Claus


ZPK = Zagreb Swimming Club ;)
___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

Re: [Scilab-users] Frequency response

[Claus Futtrup]

Hi Tim

>So, this is complicated.

I admitted from the very beginning, it's probably just me that doesn't 
know how to read the Scilab manual. It cannot be complicated for someone 
who knows and/or understand the manual.


BTW, I also notice that nobody answered the question what ZPK means? ... 
ZPK = Zagrepa Plivaci Klub ... at least the Scilab could spell out what 
it stands for.


/Claus

On 15.09.2018 22:50, Tim Wescott wrote:

So, this is complicated.

First, if you have an impulse response, in the form of a vector of
samples, then you can turn it into a FIR filter, and you can find the
frequency response of that with any of the available Scilab tools.  I'm
not sure if there's a more direct way, but if you had an impulse
response h = [1, 1/2, 1/4, 1/8], then you can turn that into H = 1 +
0.5*%z^(-1) + 0.25*%z^(-2) + 0.125*%z^(-3) (note that %z is a built-in:
%z = poly(0, 'z')).  Then you can use techniques surrounding repfreq to
get the frequency response.

Assign a sampling rate to it (i.e., by setting H.dt = 0.001, for a 1kHz
sampling rate) and you can just use the Bode plot function: bode(H).

Second, you can take your impulse response, pad it with zeros, and then
take the FFT of it -- i.e. H = fft([h zeros(1, 252)]).  You need to pad
it because the FFT is only exact for a periodic function in time --
padding it gives you an approximation of an actual impulse response
(and before someone jumps in and says you need to window it -- no, you
don't, windowing is for sample vectors of continuous signals, not
impulse responses).  But if you don't know WHY the method is
approximate, you're going to have trouble correctly massaging the data
before it's presented to the FFT, and with interpreting the results.

On Sat, 2018-09-15 at 22:15 +0200, Claus Futtrup wrote:

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(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-plo
t-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@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

Virus-free. www.avast.com
  ___
users mailing list
users@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

___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

Re: [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  on behalf of Claus 
Futtrup 

*Sent:* Saturday, September 15, 2018 11:15:08 PM
*To:* users@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@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users




<https://www.avast.com/sig-email?utm_medium=email_source=link_campaign=sig-email_content=emailclient> 
	Virus-free. www.avast.com 
<https://www.avast.com/sig-email?utm_medium=email_source=link_campaign=sig-email_content=emailclient> 





___
users mailing list
users@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
// 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 vis

Re: [Scilab-users] Frequency response

[Rafael Guerra]

Just to refresh some definitions, for non-DSP specialists like me.

For a Linear Time Invariant system with impulse response h(t) and transfer 
function H(f), input/output signals x(t) / y(t) and Fourier transforms X(f) / 
Y(f)  we have:

y(t) = x(t) * h(t)  (convolution)  ??  Y(f) = X(f).H(f)  (product)

When the input is a delta function *(t)  the output is the impulse response 
h(t) itself.

As indicated by Tim W., to obtain the system's transfer function H(f)  
("frequency response"), padding the (time) impulse response with enough zeros 
will produce enough spectral resolution in the FFT, and this is the most 
straightforward way.

Otherwise you can use for instance Scilab function's  poly, syslin, and bode;
which sounds like using a sledgehammer to crack a nut.

Regards,
Rafael

___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

Re: [Scilab-users] Frequency response

[Tim Wescott]

So, this is complicated.

First, if you have an impulse response, in the form of a vector of
samples, then you can turn it into a FIR filter, and you can find the
frequency response of that with any of the available Scilab tools.  I'm
not sure if there's a more direct way, but if you had an impulse
response h = [1, 1/2, 1/4, 1/8], then you can turn that into H = 1 +
0.5*%z^(-1) + 0.25*%z^(-2) + 0.125*%z^(-3) (note that %z is a built-in: 
%z = poly(0, 'z')).  Then you can use techniques surrounding repfreq to
get the frequency response.

Assign a sampling rate to it (i.e., by setting H.dt = 0.001, for a 1kHz
sampling rate) and you can just use the Bode plot function: bode(H).

Second, you can take your impulse response, pad it with zeros, and then
take the FFT of it -- i.e. H = fft([h zeros(1, 252)]).  You need to pad
it because the FFT is only exact for a periodic function in time --
padding it gives you an approximation of an actual impulse response
(and before someone jumps in and says you need to window it -- no, you
don't, windowing is for sample vectors of continuous signals, not
impulse responses).  But if you don't know WHY the method is
approximate, you're going to have trouble correctly massaging the data
before it's presented to the FFT, and with interpreting the results. 

On Sat, 2018-09-15 at 22:15 +0200, Claus Futtrup wrote:
> 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(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-plo
> > > t-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@lists.scilab.org
> > http://lists.scilab.org/mailman/listinfo/users
> 
>   Virus-free. www.avast.com
>  ___
> users mailing list
> users@lists.scilab.org
> http://lists.scilab.org/mailman/listinfo/users
-- 

Tim Wescott
www.wescottdesign.com
Control & Communications systems, circuit & software design.
Phone: 503.631.7815
Cell:  503.349.8432



___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

Re: [Scilab-users] Frequency response

[Rafael Guerra]

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  on behalf of Claus Futtrup 

Sent: Saturday, September 15, 2018 11:15:08 PM
To: users@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@lists.scilab.org<mailto:users@lists.scilab.org>
http://lists.scilab.org/mailman/listinfo/users



[https://ipmcdn.avast.com/images/icons/icon-envelope-tick-round-orange-animated-no-repeat-v1.gif]<https://www.avast.com/sig-email?utm_medium=email_source=link_campaign=sig-email_content=emailclient>
  Virus-free. 
www.avast.com<https://www.avast.com/sig-email?utm_medium=email_source=link_campaign=sig-email_content=emailclient>
___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

Re: [Scilab-users] Frequency response

[Claus Futtrup]

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 (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@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
___
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

Re: [Scilab-users] Frequency response

[Samuel Gougeon]

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@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users