[Scilab-users] Frequency response

classic Classic list List threaded Threaded
19 messages Options
Claus Futtrup Claus Futtrup
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] Frequency response

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);
semilogx(F,mag2db(abs(H))); grid on;
xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)');

Best regards,

Claus


Virus-free. www.avast.com

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Samuel GOUGEON Samuel GOUGEON
Reply | Threaded
Open this post in threaded view
|

Re: Frequency response

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
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Claus Futtrup Claus Futtrup
Reply | Threaded
Open this post in threaded view
|

Re: 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(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
[hidden email]
http://lists.scilab.org/mailman/listinfo/users



Virus-free. www.avast.com

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: Frequency response

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 <[hidden email]> on behalf of Claus Futtrup <[hidden email]>
Sent: Saturday, September 15, 2018 11:15:08 PM
To: [hidden email]
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(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
[hidden email]
http://lists.scilab.org/mailman/listinfo/users



Virus-free. www.avast.com

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Tim Wescott Tim Wescott
Reply | Threaded
Open this post in threaded view
|

Re: Frequency response

In reply to this post by Claus Futtrup
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
> > [hidden email]
> > http://lists.scilab.org/mailman/listinfo/users
>
> Virus-free. www.avast.com
>  _______________________________________________
> users mailing list
> [hidden email]
> 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
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: Frequency response

In reply to this post by 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 d(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
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Claus Futtrup Claus Futtrup
Reply | Threaded
Open this post in threaded view
|

Re: Frequency response

In reply to this post by Rafael Guerra
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 [hidden email] on behalf of Claus Futtrup [hidden email]
Sent: Saturday, September 15, 2018 11:15:08 PM
To: [hidden email]
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(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
[hidden email]
http://lists.scilab.org/mailman/listinfo/users



Virus-free. www.avast.com


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users



_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

stepresponse - 2018-09-16.sce (10K) Download Attachment
Claus Futtrup Claus Futtrup
Reply | Threaded
Open this post in threaded view
|

Re: Frequency response

In reply to this post by Tim Wescott
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
>>> [hidden email]
>>> http://lists.scilab.org/mailman/listinfo/users
>> Virus-free. www.avast.com
>>   _______________________________________________
>> users mailing list
>> [hidden email]
>> 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
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Alan Teeder Alan Teeder
Reply | Threaded
Open this post in threaded view
|

Re: Frequency response



-----Original Message-----
From: Claus Futtrup
Sent: Sunday, September 16, 2018 5:15 PM
To: [hidden email]
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
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Tim Wescott Tim Wescott
Reply | Threaded
Open this post in threaded view
|

Re: Frequency response

In reply to this post by Claus Futtrup
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
> > > > [hidden email]
> > > > http://lists.scilab.org/mailman/listinfo/users
> > > Virus-free. www.avast.com
> > >   _______________________________________________
> > > users mailing list
> > > [hidden email]
> > > 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
> [hidden email]
> 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
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Claus Futtrup Claus Futtrup
Reply | Threaded
Open this post in threaded view
|

Re: Frequency response

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
>>>>> [hidden email]
>>>>> http://lists.scilab.org/mailman/listinfo/users
>>>> Virus-free. www.avast.com
>>>>    _______________________________________________
>>>> users mailing list
>>>> [hidden email]
>>>> 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
>> [hidden email]
>> http://lists.scilab.org/mailman/listinfo/users


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: Frequency response

In reply to this post by Claus Futtrup

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 <[hidden email]> On Behalf Of Claus Futtrup
Sent: Sunday, September 16, 2018 5:49 PM
To: [hidden email]
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 [hidden email] on behalf of Claus Futtrup [hidden email]
Sent: Saturday, September 15, 2018 11:15:08 PM
To: [hidden email]
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(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
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

 

 

Virus-free. www.avast.com




_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users

 


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Christophe Dang Ngoc Chan Christophe Dang Ngoc Chan
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} Re: Frequency response

In reply to this post by Claus Futtrup
Hello,

> De : users [mailto:[hidden email]] De la part de Claus Futtrup
> Envoyé : dimanche 16 septembre 2018 19:14
>
> 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.

I think you got it right.

My qwant-fu with "zpk zero pole" shows up as first answer
mathworks.com/help/control/ref/zpk.html

"continuous-time zero-pole-gain model with zeros Z, poles P, and gain(s) K."

Regards

--
Christophe Dang Ngoc Chan
Mechanical calculation engineer
This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Tan Chin Luh-2 Tan Chin Luh-2
Reply | Threaded
Open this post in threaded view
|

Re: Frequency response

In reply to this post by Rafael Guerra
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
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: Frequency response

Chin Luh:  good to know but how does that solve Claus Futtrup specific problem?
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} Re: Frequency response

In reply to this post by Christophe Dang Ngoc Chan
Redundant qwant-fu'ing ...
response was already provided


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Tan Chin Luh-2 Tan Chin Luh-2
Reply | Threaded
Open this post in threaded view
|

Re: Frequency response

In reply to this post by Rafael Guerra
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
[hidden email]
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
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: Frequency response

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
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Claus Futtrup Claus Futtrup
Reply | Threaded
Open this post in threaded view
|

Re: Frequency response

In reply to this post by Tan Chin Luh-2
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
[hidden email]
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
[hidden email]
http://lists.scilab.org/mailman/listinfo/users



Virus-free. www.avast.com

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users