# [Scilab-users] Frequency response

20 messages
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?

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
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
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?

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
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?

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
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
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
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?

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

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
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
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
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
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])`
` `

`[ [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?

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
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
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
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
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
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
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
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:

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