Menu
▾
▴
octave-cvsupdate
|
From: <ad...@us...> - 2008-07-30 08:47:40
|
Revision: 5206
http://octave.svn.sourceforge.net/octave/?rev=5206&view=rev
Author: adb014
Date: 2008-07-30 08:47:50 +0000 (Wed, 30 Jul 2008)
Log Message:
-----------
Add 1 to phase in upsample/downsample for matlab compatiblity (For Jonathan Coker)
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/downsample.m
trunk/octave-forge/main/signal/inst/upsample.m
Modified: trunk/octave-forge/main/signal/inst/downsample.m
===================================================================
--- trunk/octave-forge/main/signal/inst/downsample.m 2008-07-30 07:59:40 UTC (rev 5205)
+++ trunk/octave-forge/main/signal/inst/downsample.m 2008-07-30 08:47:50 UTC (rev 5206)
@@ -7,8 +7,8 @@
## it prefilters the high frequency components of the signal and
## avoids aliasing effects.
##
-## @deftypefnx {Function File} @var{y} = downsample(@var{x},@var{n},@var{phase})
-## Select every nth element starting at sample @var{phase}.
+## @deftypefnx {Function File} @var{y} = downsample(@var{x},@var{n},@var{offset})
+## Select every nth element starting at sample @var{offset}.
## @end deftypefn
## @seealso{decimate, interp, resample, upfirdn, upsample}
@@ -16,24 +16,24 @@
## This function is public domain
function y = downsample(x,n,phase)
- if nargin<2 || nargin>3, usage('downsample(x,n,[phase]'); end
- if nargin==2, phase = 1; end
+ if nargin<2 || nargin>3, usage('downsample(x,n,[offset]'); end
+ if nargin==2, phase = 0; end
- if phase > n
+ if phase > n - 1
warning("This is incompatible with Matlab (phase = 0:n-1). See\
octave-forge signal package release notes for details." )
end
if isvector(x)
- y = x(phase:n:end);
+ y = x(phase + 1:n:end);
else
- y = x(phase:n:end,:);
+ y = x(phase + 1:n:end,:);
end
end
%!assert(downsample([1,2,3,4,5],2),[1,3,5]);
%!assert(downsample([1;2;3;4;5],2),[1;3;5]);
%!assert(downsample([1,2;3,4;5,6;7,8;9,10],2),[1,2;5,6;9,10]);
-%!assert(downsample([1,2,3,4,5],2,2),[2,4]);
-%!assert(downsample([1,2;3,4;5,6;7,8;9,10],2,2),[3,4;7,8]);
+%!assert(downsample([1,2,3,4,5],2,1),[2,4]);
+%!assert(downsample([1,2;3,4;5,6;7,8;9,10],2,1),[3,4;7,8]);
Modified: trunk/octave-forge/main/signal/inst/upsample.m
===================================================================
--- trunk/octave-forge/main/signal/inst/upsample.m 2008-07-30 07:59:40 UTC (rev 5205)
+++ trunk/octave-forge/main/signal/inst/upsample.m 2008-07-30 08:47:50 UTC (rev 5206)
@@ -3,7 +3,7 @@
## Upsample the signal, inserting n-1 zeros between every element.
## If @var{x} is a matrix, upsample every column.
##
-## @deftypefnx {Function File} @var{y} = upsample(@var{x},@var{n},@var{phase})
+## @deftypefnx {Function File} @var{y} = upsample(@var{x},@var{n},@var{offset})
## Control the position of the inserted sample in the block of n zeros.
## @end deftypefn
## @seealso{decimate, downsample, interp, resample, upfirdn}
@@ -13,22 +13,27 @@
function y = upsample(x,n,phase)
if nargin<2 || nargin>3, usage('upsample(x,n,[phase]'); end
- if nargin==2, phase = 1; end
+ if nargin==2, phase = 0; end
+ if phase > n - 1
+ warning("This is incompatible with Matlab (phase = 0:n-1). See\
+ octave-forge signal package release notes for details." )
+ end
+
[nr,nc] = size(x);
if any([nr,nc]==1),
y = zeros(n*nr*nc,1);
- y(phase:n:end) = x;
+ y(phase + 1:n:end) = x;
if nr==1, y = y.'; end
else
y = zeros(n*nr,nc);
- y(phase:n:end,:) = x;
+ y(phase + 1:n:end,:) = x;
end
end
%!assert(upsample([1,3,5],2),[1,0,3,0,5,0]);
%!assert(upsample([1;3;5],2),[1;0;3;0;5;0]);
%!assert(upsample([1,2;5,6;9,10],2),[1,2;0,0;5,6;0,0;9,10;0,0]);
-%!assert(upsample([2,4],2,2),[0,2,0,4]);
-%!assert(upsample([3,4;7,8],2,2),[0,0;3,4;0,0;7,8]);
+%!assert(upsample([2,4],2,1),[0,2,0,4]);
+%!assert(upsample([3,4;7,8],2,1),[0,0;3,4;0,0;7,8]);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <pki...@us...> - 2008-10-21 23:22:06
|
Revision: 5384
http://octave.svn.sourceforge.net/octave/?rev=5384&view=rev
Author: pkienzle
Date: 2008-10-21 23:22:00 +0000 (Tue, 21 Oct 2008)
Log Message:
-----------
Use Fs rather than T for sampling frequency in IIR filter design functions
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/bilinear.m
trunk/octave-forge/main/signal/inst/butter.m
trunk/octave-forge/main/signal/inst/buttord.m
trunk/octave-forge/main/signal/inst/cheb1ord.m
trunk/octave-forge/main/signal/inst/cheb2ord.m
trunk/octave-forge/main/signal/inst/cheby1.m
trunk/octave-forge/main/signal/inst/cheby2.m
trunk/octave-forge/main/signal/inst/ellip.m
trunk/octave-forge/main/signal/inst/ellipord.m
Modified: trunk/octave-forge/main/signal/inst/bilinear.m
===================================================================
--- trunk/octave-forge/main/signal/inst/bilinear.m 2008-10-20 17:47:58 UTC (rev 5383)
+++ trunk/octave-forge/main/signal/inst/bilinear.m 2008-10-21 23:22:00 UTC (rev 5384)
@@ -13,20 +13,20 @@
## You should have received a copy of the GNU General Public License
## along with this program; If not, see <http://www.gnu.org/licenses/>.
-## usage: [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T)
-## [Zb, Za] = bilinear(Sb, Sa, T)
+## usage: [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, Fs)
+## [Zb, Za] = bilinear(Sb, Sa, Fs)
##
## Transform a s-plane filter specification into a z-plane
## specification. Filters can be specified in either zero-pole-gain or
## transfer function form. The input form does not have to match the
-## output form. T is the sampling frequency represented in the z plane.
+## output form. Fs is the sampling frequency represented in the z plane.
##
## Theory: Given a piecewise flat filter design, you can transform it
## from the s-plane to the z-plane while maintaining the band edges by
## means of the bilinear transform. This maps the left hand side of the
## s-plane into the interior of the unit circle. The mapping is highly
## non-linear, so you must design your filter with band edges in the
-## s-plane positioned at 2/T tan(w*T/2) so that they will be positioned
+## s-plane positioned at 2/Fs tan(w*Fs/2) so that they will be positioned
## at w after the bilinear transform is complete.
##
## The following table summarizes the transformation:
@@ -35,9 +35,9 @@
## | Transform | Zero at x | Pole at x |
## | H(S) | H(S) = S-x | H(S)=1/(S-x) |
## +---------------+-----------------------+----------------------+
-## | 2 z-1 | zero: (2+xT)/(2-xT) | zero: -1 |
-## | S -> - --- | pole: -1 | pole: (2+xT)/(2-xT) |
-## | T z+1 | gain: (2-xT)/T | gain: (2-xT)/T |
+## | 2 z-1 | zero: (2+xFs)/(2-xFs) | zero: -1 |
+## | S -> -- --- | pole: -1 | pole: (2+xFs)/(2-xFs)|
+## | Fs z+1 | gain: (2-xFs)/Fs | gain: (2-xFs)/Fs |
## +---------------+-----------------------+----------------------+
##
## With tedious algebra, you can derive the above formulae yourself by
@@ -69,13 +69,13 @@
## Author: Paul Kienzle <pki...@us...>
-function [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T)
+function [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, Fs)
if nargin==3
- T = Sg;
+ Fs = Sg;
[Sz, Sp, Sg] = tf2zp(Sz, Sp);
elseif nargin!=4
- usage("[Zz, Zp, Zg]=bilinear(Sz,Sp,Sg,T) or [Zb, Za]=blinear(Sb,Sa,T)");
+ usage("[Zz, Zp, Zg]=bilinear(Sz,Sp,Sg,Fs) or [Zb, Za]=blinear(Sb,Sa,Fs)");
end;
p = length(Sp);
@@ -85,17 +85,17 @@
end
## ---------------- ------------------------- ------------------------
-## Bilinear zero: (2+xT)/(2-xT) pole: (2+xT)/(2-xT)
-## 2 z-1 pole: -1 zero: -1
-## S -> - --- gain: (2-xT)/T gain: (2-xT)/T
-## T z+1
+## Bilinear zero: (2+xFs)/(2-xFs) pole: (2+xFs)/(2-xFs)
+## 2 z-1 pole: -1 zero: -1
+## S -> -- --- gain: (2-xFs)/Fs gain: (2-xFs)/Fs
+## Fs z+1
## ---------------- ------------------------- ------------------------
- Zg = real(Sg * prod((2-Sz*T)/T) / prod((2-Sp*T)/T));
- Zp = (2+Sp*T)./(2-Sp*T);
+ Zg = real(Sg * prod((2-Sz*Fs)/Fs) / prod((2-Sp*Fs)/Fs));
+ Zp = (2+Sp*Fs)./(2-Sp*Fs);
if isempty(Sz)
Zz = -ones(size(Zp));
else
- Zz = [(2+Sz*T)./(2-Sz*T)];
+ Zz = [(2+Sz*Fs)./(2-Sz*Fs)];
Zz = postpad(Zz, p, -1);
end
Modified: trunk/octave-forge/main/signal/inst/butter.m
===================================================================
--- trunk/octave-forge/main/signal/inst/butter.m 2008-10-20 17:47:58 UTC (rev 5383)
+++ trunk/octave-forge/main/signal/inst/butter.m 2008-10-21 23:22:00 UTC (rev 5384)
@@ -87,8 +87,8 @@
## Prewarp to the band edges to s plane
if digital
- T = 2; # sampling frequency of 2 Hz
- W = 2/T*tan(pi*W/T);
+ Fs = 2; # sampling frequency of 2 Hz
+ W = 2/Fs*tan(pi*W/Fs);
endif
## Generate splane poles for the prototype butterworth filter
@@ -104,7 +104,7 @@
## Use bilinear transform to convert poles to the z plane
if digital
- [zero, pole, gain] = bilinear(zero, pole, gain, T);
+ [zero, pole, gain] = bilinear(zero, pole, gain, Fs);
endif
## convert to the correct output form
Modified: trunk/octave-forge/main/signal/inst/buttord.m
===================================================================
--- trunk/octave-forge/main/signal/inst/buttord.m 2008-10-20 17:47:58 UTC (rev 5383)
+++ trunk/octave-forge/main/signal/inst/buttord.m 2008-10-21 23:22:00 UTC (rev 5384)
@@ -55,7 +55,7 @@
warning("buttord: seems to overdesign bandpass and bandreject filters");
end
- T = 2;
+ Fs = 2;
## if high pass, reverse the sense of the test
stop = find(Wp > Ws);
@@ -63,8 +63,8 @@
Ws(stop) = 1-Ws(stop); # subtract from ones(1,length(stop))
## warp the target frequencies according to the bilinear transform
- Ws = (2/T)*tan(pi*Ws./T);
- Wp = (2/T)*tan(pi*Wp./T);
+ Ws = (2/Fs)*tan(pi*Ws./Fs);
+ Wp = (2/Fs)*tan(pi*Wp./Fs);
## compute minimum n which satisfies all band edge conditions
## the factor 1/length(Wp) is an artificial correction for the
@@ -77,7 +77,7 @@
Wc = exp(log(Wp) - qp/2/n);
## unwarp the returned frequency
- Wc = atan(T/2*Wc)*T/pi;
+ Wc = atan(Fs/2*Wc)*Fs/pi;
## if high pass, reverse the sense of the test
Wc(stop) = 1-Wc(stop);
Modified: trunk/octave-forge/main/signal/inst/cheb1ord.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cheb1ord.m 2008-10-20 17:47:58 UTC (rev 5383)
+++ trunk/octave-forge/main/signal/inst/cheb1ord.m 2008-10-21 23:22:00 UTC (rev 5384)
@@ -44,14 +44,14 @@
error("cheb1ord: Wp(1)<Ws(1)<Ws(2)<Wp(2) or Ws(1)<Wp(1)<Wp(2)<Ws(2)");
end
- T = 2;
+ Fs = 2;
## returned frequency is the same as the input frequency
Wc = Wp;
## warp the target frequencies according to the bilinear transform
- Ws = (2/T)*tan(pi*Ws./T);
- Wp = (2/T)*tan(pi*Wp./T);
+ Ws = (2/Fs)*tan(pi*Ws./Fs);
+ Wp = (2/Fs)*tan(pi*Wp./Fs);
if (Wp(1) < Ws(1))
## low pass
Modified: trunk/octave-forge/main/signal/inst/cheb2ord.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cheb2ord.m 2008-10-20 17:47:58 UTC (rev 5383)
+++ trunk/octave-forge/main/signal/inst/cheb2ord.m 2008-10-21 23:22:00 UTC (rev 5384)
@@ -44,14 +44,14 @@
error("cheb2ord: Wp(1)<Ws(1)<Ws(2)<Wp(2) or Ws(1)<Wp(1)<Wp(2)<Ws(2)");
end
- T = 2;
+ Fs = 2;
## returned frequency is the same as the input frequency
Wc = Ws;
## warp the target frequencies according to the bilinear transform
- Ws = (2/T)*tan(pi*Ws./T);
- Wp = (2/T)*tan(pi*Wp./T);
+ Ws = (2/Fs)*tan(pi*Ws./Fs);
+ Wp = (2/Fs)*tan(pi*Wp./Fs);
if (Wp(1) < Ws(1))
## low pass
Modified: trunk/octave-forge/main/signal/inst/cheby1.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cheby1.m 2008-10-20 17:47:58 UTC (rev 5383)
+++ trunk/octave-forge/main/signal/inst/cheby1.m 2008-10-21 23:22:00 UTC (rev 5384)
@@ -89,8 +89,8 @@
## Prewarp to the band edges to s plane
if digital
- T = 2; # sampling frequency of 2 Hz
- W = 2/T*tan(pi*W/T);
+ Fs = 2; # sampling frequency of 2 Hz
+ W = 2/Fs*tan(pi*W/Fs);
endif
## Generate splane poles and zeros for the chebyshev type 1 filter
@@ -114,7 +114,7 @@
## Use bilinear transform to convert poles to the z plane
if digital
- [zero, pole, gain] = bilinear(zero, pole, gain, T);
+ [zero, pole, gain] = bilinear(zero, pole, gain, Fs);
endif
## convert to the correct output form
Modified: trunk/octave-forge/main/signal/inst/cheby2.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cheby2.m 2008-10-20 17:47:58 UTC (rev 5383)
+++ trunk/octave-forge/main/signal/inst/cheby2.m 2008-10-21 23:22:00 UTC (rev 5384)
@@ -13,7 +13,7 @@
## You should have received a copy of the GNU General Public License
## along with this program; If not, see <http://www.gnu.org/licenses/>.
-## Generate an Chebyshev type II filter with Rs dB of stop band attenuation.
+## Generate an Chebyshev type II filter with Rs dB of stop band ripple.
##
## [b, a] = cheby2(n, Rs, Wc)
## low pass filter with cutoff pi*Wc radians
@@ -90,8 +90,8 @@
## Prewarp to the band edges to s plane
if digital
- T = 2; # sampling frequency of 2 Hz
- W = 2/T*tan(pi*W/T);
+ Fs = 2; # sampling frequency of 2 Hz
+ W = 2/Fs*tan(pi*W/Fs);
endif
## Generate splane poles and zeros for the chebyshev type 2 filter
@@ -123,7 +123,7 @@
## Use bilinear transform to convert poles to the z plane
if digital
- [zero, pole, gain] = bilinear(zero, pole, gain, T);
+ [zero, pole, gain] = bilinear(zero, pole, gain, Fs);
endif
## convert to the correct output form
Modified: trunk/octave-forge/main/signal/inst/ellip.m
===================================================================
--- trunk/octave-forge/main/signal/inst/ellip.m 2008-10-20 17:47:58 UTC (rev 5383)
+++ trunk/octave-forge/main/signal/inst/ellip.m 2008-10-21 23:22:00 UTC (rev 5384)
@@ -102,8 +102,8 @@
##Prewarp the digital frequencies
if digital
- T = 2; # sampling frequency of 2 Hz
- W = tan(pi*W/T);
+ Fs = 2; # sampling frequency of 2 Hz
+ W = tan(pi*W/Fs);
endif
##Generate s-plane poles, zeros and gain
@@ -114,7 +114,7 @@
## Use bilinear transform to convert poles to the z plane
if digital
- [zero, pole, gain] = bilinear(zero, pole, gain, T);
+ [zero, pole, gain] = bilinear(zero, pole, gain, Fs);
endif
Modified: trunk/octave-forge/main/signal/inst/ellipord.m
===================================================================
--- trunk/octave-forge/main/signal/inst/ellipord.m 2008-10-20 17:47:58 UTC (rev 5383)
+++ trunk/octave-forge/main/signal/inst/ellipord.m 2008-10-21 23:22:00 UTC (rev 5384)
@@ -55,10 +55,10 @@
# sampling frequency of 2 Hz
-T = 2;
+Fs = 2;
-Wpw = tan(pi.*Wp./T); # prewarp
-Wsw = tan(pi.*Ws./T); # prewarp
+Wpw = tan(pi.*Wp./Fs); # prewarp
+Wsw = tan(pi.*Ws./Fs); # prewarp
##pass/stop band to low pass filter transform:
if (length(Wpw)==2 && length(Wsw)==2)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <pki...@us...> - 2008-10-27 10:01:15
|
Revision: 5386
http://octave.svn.sourceforge.net/octave/?rev=5386&view=rev
Author: pkienzle
Date: 2008-10-27 09:56:54 +0000 (Mon, 27 Oct 2008)
Log Message:
-----------
signal: undoing change from T->Fs in iir filters
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/bilinear.m
trunk/octave-forge/main/signal/inst/butter.m
trunk/octave-forge/main/signal/inst/buttord.m
trunk/octave-forge/main/signal/inst/cheb1ord.m
trunk/octave-forge/main/signal/inst/cheb2ord.m
trunk/octave-forge/main/signal/inst/cheby1.m
trunk/octave-forge/main/signal/inst/cheby2.m
trunk/octave-forge/main/signal/inst/ellip.m
trunk/octave-forge/main/signal/inst/ellipord.m
Modified: trunk/octave-forge/main/signal/inst/bilinear.m
===================================================================
--- trunk/octave-forge/main/signal/inst/bilinear.m 2008-10-25 09:09:50 UTC (rev 5385)
+++ trunk/octave-forge/main/signal/inst/bilinear.m 2008-10-27 09:56:54 UTC (rev 5386)
@@ -13,20 +13,20 @@
## You should have received a copy of the GNU General Public License
## along with this program; If not, see <http://www.gnu.org/licenses/>.
-## usage: [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, Fs)
-## [Zb, Za] = bilinear(Sb, Sa, Fs)
+## usage: [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T)
+## [Zb, Za] = bilinear(Sb, Sa, T)
##
## Transform a s-plane filter specification into a z-plane
## specification. Filters can be specified in either zero-pole-gain or
## transfer function form. The input form does not have to match the
-## output form. Fs is the sampling frequency represented in the z plane.
+## output form. T is the sampling frequency represented in the z plane.
##
## Theory: Given a piecewise flat filter design, you can transform it
## from the s-plane to the z-plane while maintaining the band edges by
## means of the bilinear transform. This maps the left hand side of the
## s-plane into the interior of the unit circle. The mapping is highly
## non-linear, so you must design your filter with band edges in the
-## s-plane positioned at 2/Fs tan(w*Fs/2) so that they will be positioned
+## s-plane positioned at 2/T tan(w*T/2) so that they will be positioned
## at w after the bilinear transform is complete.
##
## The following table summarizes the transformation:
@@ -35,9 +35,9 @@
## | Transform | Zero at x | Pole at x |
## | H(S) | H(S) = S-x | H(S)=1/(S-x) |
## +---------------+-----------------------+----------------------+
-## | 2 z-1 | zero: (2+xFs)/(2-xFs) | zero: -1 |
-## | S -> -- --- | pole: -1 | pole: (2+xFs)/(2-xFs)|
-## | Fs z+1 | gain: (2-xFs)/Fs | gain: (2-xFs)/Fs |
+## | 2 z-1 | zero: (2+xT)/(2-xT) | zero: -1 |
+## | S -> - --- | pole: -1 | pole: (2+xT)/(2-xT) |
+## | T z+1 | gain: (2-xT)/T | gain: (2-xT)/T |
## +---------------+-----------------------+----------------------+
##
## With tedious algebra, you can derive the above formulae yourself by
@@ -69,13 +69,13 @@
## Author: Paul Kienzle <pki...@us...>
-function [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, Fs)
+function [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T)
if nargin==3
- Fs = Sg;
+ T = Sg;
[Sz, Sp, Sg] = tf2zp(Sz, Sp);
elseif nargin!=4
- usage("[Zz, Zp, Zg]=bilinear(Sz,Sp,Sg,Fs) or [Zb, Za]=blinear(Sb,Sa,Fs)");
+ usage("[Zz, Zp, Zg]=bilinear(Sz,Sp,Sg,T) or [Zb, Za]=blinear(Sb,Sa,T)");
end;
p = length(Sp);
@@ -85,17 +85,17 @@
end
## ---------------- ------------------------- ------------------------
-## Bilinear zero: (2+xFs)/(2-xFs) pole: (2+xFs)/(2-xFs)
-## 2 z-1 pole: -1 zero: -1
-## S -> -- --- gain: (2-xFs)/Fs gain: (2-xFs)/Fs
-## Fs z+1
+## Bilinear zero: (2+xT)/(2-xT) pole: (2+xT)/(2-xT)
+## 2 z-1 pole: -1 zero: -1
+## S -> - --- gain: (2-xT)/T gain: (2-xT)/T
+## T z+1
## ---------------- ------------------------- ------------------------
- Zg = real(Sg * prod((2-Sz*Fs)/Fs) / prod((2-Sp*Fs)/Fs));
- Zp = (2+Sp*Fs)./(2-Sp*Fs);
+ Zg = real(Sg * prod((2-Sz*T)/T) / prod((2-Sp*T)/T));
+ Zp = (2+Sp*T)./(2-Sp*T);
if isempty(Sz)
Zz = -ones(size(Zp));
else
- Zz = [(2+Sz*Fs)./(2-Sz*Fs)];
+ Zz = [(2+Sz*T)./(2-Sz*T)];
Zz = postpad(Zz, p, -1);
end
Modified: trunk/octave-forge/main/signal/inst/butter.m
===================================================================
--- trunk/octave-forge/main/signal/inst/butter.m 2008-10-25 09:09:50 UTC (rev 5385)
+++ trunk/octave-forge/main/signal/inst/butter.m 2008-10-27 09:56:54 UTC (rev 5386)
@@ -87,8 +87,8 @@
## Prewarp to the band edges to s plane
if digital
- Fs = 2; # sampling frequency of 2 Hz
- W = 2/Fs*tan(pi*W/Fs);
+ T = 2; # sampling frequency of 2 Hz
+ W = 2/T*tan(pi*W/T);
endif
## Generate splane poles for the prototype butterworth filter
@@ -104,7 +104,7 @@
## Use bilinear transform to convert poles to the z plane
if digital
- [zero, pole, gain] = bilinear(zero, pole, gain, Fs);
+ [zero, pole, gain] = bilinear(zero, pole, gain, T);
endif
## convert to the correct output form
Modified: trunk/octave-forge/main/signal/inst/buttord.m
===================================================================
--- trunk/octave-forge/main/signal/inst/buttord.m 2008-10-25 09:09:50 UTC (rev 5385)
+++ trunk/octave-forge/main/signal/inst/buttord.m 2008-10-27 09:56:54 UTC (rev 5386)
@@ -55,7 +55,7 @@
warning("buttord: seems to overdesign bandpass and bandreject filters");
end
- Fs = 2;
+ T = 2;
## if high pass, reverse the sense of the test
stop = find(Wp > Ws);
@@ -63,8 +63,8 @@
Ws(stop) = 1-Ws(stop); # subtract from ones(1,length(stop))
## warp the target frequencies according to the bilinear transform
- Ws = (2/Fs)*tan(pi*Ws./Fs);
- Wp = (2/Fs)*tan(pi*Wp./Fs);
+ Ws = (2/T)*tan(pi*Ws./T);
+ Wp = (2/T)*tan(pi*Wp./T);
## compute minimum n which satisfies all band edge conditions
## the factor 1/length(Wp) is an artificial correction for the
@@ -77,7 +77,7 @@
Wc = exp(log(Wp) - qp/2/n);
## unwarp the returned frequency
- Wc = atan(Fs/2*Wc)*Fs/pi;
+ Wc = atan(T/2*Wc)*T/pi;
## if high pass, reverse the sense of the test
Wc(stop) = 1-Wc(stop);
Modified: trunk/octave-forge/main/signal/inst/cheb1ord.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cheb1ord.m 2008-10-25 09:09:50 UTC (rev 5385)
+++ trunk/octave-forge/main/signal/inst/cheb1ord.m 2008-10-27 09:56:54 UTC (rev 5386)
@@ -44,14 +44,14 @@
error("cheb1ord: Wp(1)<Ws(1)<Ws(2)<Wp(2) or Ws(1)<Wp(1)<Wp(2)<Ws(2)");
end
- Fs = 2;
+ T = 2;
## returned frequency is the same as the input frequency
Wc = Wp;
## warp the target frequencies according to the bilinear transform
- Ws = (2/Fs)*tan(pi*Ws./Fs);
- Wp = (2/Fs)*tan(pi*Wp./Fs);
+ Ws = (2/T)*tan(pi*Ws./T);
+ Wp = (2/T)*tan(pi*Wp./T);
if (Wp(1) < Ws(1))
## low pass
Modified: trunk/octave-forge/main/signal/inst/cheb2ord.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cheb2ord.m 2008-10-25 09:09:50 UTC (rev 5385)
+++ trunk/octave-forge/main/signal/inst/cheb2ord.m 2008-10-27 09:56:54 UTC (rev 5386)
@@ -44,14 +44,14 @@
error("cheb2ord: Wp(1)<Ws(1)<Ws(2)<Wp(2) or Ws(1)<Wp(1)<Wp(2)<Ws(2)");
end
- Fs = 2;
+ T = 2;
## returned frequency is the same as the input frequency
Wc = Ws;
## warp the target frequencies according to the bilinear transform
- Ws = (2/Fs)*tan(pi*Ws./Fs);
- Wp = (2/Fs)*tan(pi*Wp./Fs);
+ Ws = (2/T)*tan(pi*Ws./T);
+ Wp = (2/T)*tan(pi*Wp./T);
if (Wp(1) < Ws(1))
## low pass
Modified: trunk/octave-forge/main/signal/inst/cheby1.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cheby1.m 2008-10-25 09:09:50 UTC (rev 5385)
+++ trunk/octave-forge/main/signal/inst/cheby1.m 2008-10-27 09:56:54 UTC (rev 5386)
@@ -89,8 +89,8 @@
## Prewarp to the band edges to s plane
if digital
- Fs = 2; # sampling frequency of 2 Hz
- W = 2/Fs*tan(pi*W/Fs);
+ T = 2; # sampling frequency of 2 Hz
+ W = 2/T*tan(pi*W/T);
endif
## Generate splane poles and zeros for the chebyshev type 1 filter
@@ -114,7 +114,7 @@
## Use bilinear transform to convert poles to the z plane
if digital
- [zero, pole, gain] = bilinear(zero, pole, gain, Fs);
+ [zero, pole, gain] = bilinear(zero, pole, gain, T);
endif
## convert to the correct output form
Modified: trunk/octave-forge/main/signal/inst/cheby2.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cheby2.m 2008-10-25 09:09:50 UTC (rev 5385)
+++ trunk/octave-forge/main/signal/inst/cheby2.m 2008-10-27 09:56:54 UTC (rev 5386)
@@ -13,7 +13,7 @@
## You should have received a copy of the GNU General Public License
## along with this program; If not, see <http://www.gnu.org/licenses/>.
-## Generate an Chebyshev type II filter with Rs dB of stop band ripple.
+## Generate an Chebyshev type II filter with Rs dB of stop band attenuation.
##
## [b, a] = cheby2(n, Rs, Wc)
## low pass filter with cutoff pi*Wc radians
@@ -90,8 +90,8 @@
## Prewarp to the band edges to s plane
if digital
- Fs = 2; # sampling frequency of 2 Hz
- W = 2/Fs*tan(pi*W/Fs);
+ T = 2; # sampling frequency of 2 Hz
+ W = 2/T*tan(pi*W/T);
endif
## Generate splane poles and zeros for the chebyshev type 2 filter
@@ -123,7 +123,7 @@
## Use bilinear transform to convert poles to the z plane
if digital
- [zero, pole, gain] = bilinear(zero, pole, gain, Fs);
+ [zero, pole, gain] = bilinear(zero, pole, gain, T);
endif
## convert to the correct output form
Modified: trunk/octave-forge/main/signal/inst/ellip.m
===================================================================
--- trunk/octave-forge/main/signal/inst/ellip.m 2008-10-25 09:09:50 UTC (rev 5385)
+++ trunk/octave-forge/main/signal/inst/ellip.m 2008-10-27 09:56:54 UTC (rev 5386)
@@ -102,8 +102,8 @@
##Prewarp the digital frequencies
if digital
- Fs = 2; # sampling frequency of 2 Hz
- W = tan(pi*W/Fs);
+ T = 2; # sampling frequency of 2 Hz
+ W = tan(pi*W/T);
endif
##Generate s-plane poles, zeros and gain
@@ -114,7 +114,7 @@
## Use bilinear transform to convert poles to the z plane
if digital
- [zero, pole, gain] = bilinear(zero, pole, gain, Fs);
+ [zero, pole, gain] = bilinear(zero, pole, gain, T);
endif
Modified: trunk/octave-forge/main/signal/inst/ellipord.m
===================================================================
--- trunk/octave-forge/main/signal/inst/ellipord.m 2008-10-25 09:09:50 UTC (rev 5385)
+++ trunk/octave-forge/main/signal/inst/ellipord.m 2008-10-27 09:56:54 UTC (rev 5386)
@@ -55,10 +55,10 @@
# sampling frequency of 2 Hz
-Fs = 2;
+T = 2;
-Wpw = tan(pi.*Wp./Fs); # prewarp
-Wsw = tan(pi.*Ws./Fs); # prewarp
+Wpw = tan(pi.*Wp./T); # prewarp
+Wsw = tan(pi.*Ws./T); # prewarp
##pass/stop band to low pass filter transform:
if (length(Wpw)==2 && length(Wsw)==2)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ha...@us...> - 2009-06-14 18:06:29
|
Revision: 5940
http://octave.svn.sourceforge.net/octave/?rev=5940&view=rev
Author: hauberg
Date: 2009-06-14 18:06:28 +0000 (Sun, 14 Jun 2009)
Log Message:
-----------
Sign error in documentation (from Torsten Finke)
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/dct.m
trunk/octave-forge/main/signal/inst/idct.m
Modified: trunk/octave-forge/main/signal/inst/dct.m
===================================================================
--- trunk/octave-forge/main/signal/inst/dct.m 2009-06-13 22:41:22 UTC (rev 5939)
+++ trunk/octave-forge/main/signal/inst/dct.m 2009-06-14 18:06:28 UTC (rev 5940)
@@ -23,7 +23,7 @@
## The discrete cosine transform X of x can be defined as follows:
##
## N-1
-## X[k] = w(k) sum x[n] cos (pi (2n-1) k / 2N ), k = 0, ..., N-1
+## X[k] = w(k) sum x[n] cos (pi (2n+1) k / 2N ), k = 0, ..., N-1
## n=0
##
## with w(0) = sqrt(1/N) and w(k) = sqrt(2/N), k = 1, ..., N-1. There
@@ -37,7 +37,7 @@
## the discrete cosine transform of x at k is as follows:
##
## N-1
-## X[k] = sum 2 x[n] cos (pi (2n-1) k / 2N )
+## X[k] = sum 2 x[n] cos (pi (2n+1) k / 2N )
## n=0
##
## which can be computed using:
Modified: trunk/octave-forge/main/signal/inst/idct.m
===================================================================
--- trunk/octave-forge/main/signal/inst/idct.m 2009-06-13 22:41:22 UTC (rev 5939)
+++ trunk/octave-forge/main/signal/inst/idct.m 2009-06-14 18:06:28 UTC (rev 5940)
@@ -23,7 +23,7 @@
## The inverse discrete cosine transform x of X can be defined as follows:
##
## N-1
-## x[n] = sum w(k) X[k] cos (pi (2n-1) k / 2N ), k = 0, ..., N-1
+## x[n] = sum w(k) X[k] cos (pi (2n+1) k / 2N ), k = 0, ..., N-1
## k=0
##
## with w(0) = sqrt(1/N) and w(k) = sqrt(2/N), k = 1, ..., N-1
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ha...@us...> - 2009-07-18 21:06:20
|
Revision: 6026
http://octave.svn.sourceforge.net/octave/?rev=6026&view=rev
Author: hauberg
Date: 2009-07-18 21:06:08 +0000 (Sat, 18 Jul 2009)
Log Message:
-----------
New functions: besself and besselap (from Thomas Sailer)
Added Paths:
-----------
trunk/octave-forge/main/signal/inst/besselap.m
trunk/octave-forge/main/signal/inst/besself.m
Added: trunk/octave-forge/main/signal/inst/besselap.m
===================================================================
--- trunk/octave-forge/main/signal/inst/besselap.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/besselap.m 2009-07-18 21:06:08 UTC (rev 6026)
@@ -0,0 +1,52 @@
+## Copyright (C) 2009 Thomas Sailer
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+## Return bessel analog filter prototype.
+##
+## References:
+##
+## http://en.wikipedia.org/wiki/Bessel_polynomials
+
+function [zero, pole, gain] = besselap (n)
+
+ if (nargin>1 || nargin<1)
+ usage ("[z, p, g] = besselap (n)");
+ end
+
+ ## interpret the input parameters
+ if (!(length(n)==1 && n == round(n) && n > 0))
+ error ("besselap: filter order n must be a positive integer");
+ end
+
+ p0=1;
+ p1=[1 1];
+ for nn=2:n;
+ px=(2*nn-1)*p1;
+ py=[p0 0 0];
+ px=prepad(px,max(length(px),length(py)),0);
+ py=prepad(py,length(px));
+ p0=p1;
+ p1=px+py;
+ endfor;
+ % p1 now contains the reverse bessel polynomial for n
+
+ % scale it by replacing s->s/w0 so that the gain becomes 1
+ p1=p1.*p1(length(p1)).^((length(p1)-1:-1:0)/(length(p1)-1));
+
+ zero=[];
+ pole=roots(p1);
+ gain=1;
+
+endfunction
Added: trunk/octave-forge/main/signal/inst/besself.m
===================================================================
--- trunk/octave-forge/main/signal/inst/besself.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/besself.m 2009-07-18 21:06:08 UTC (rev 6026)
@@ -0,0 +1,119 @@
+## Copyright (C) 2009 Thomas Sailer
+## Copyright (C) 1999 Paul Kienzle
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+## Generate a bessel filter.
+## Default is a Laplace space (s) filter.
+##
+## [b,a] = besself(n, Wc)
+## low pass filter with cutoff pi*Wc radians
+##
+## [b,a] = besself(n, Wc, 'high')
+## high pass filter with cutoff pi*Wc radians
+##
+## [b,a] = besself(n, [Wl, Wh])
+## band pass filter with edges pi*Wl and pi*Wh radians
+##
+## [b,a] = besself(n, [Wl, Wh], 'stop')
+## band reject filter with edges pi*Wl and pi*Wh radians
+##
+## [z,p,g] = besself(...)
+## return filter as zero-pole-gain rather than coefficients of the
+## numerator and denominator polynomials.
+##
+## [...] = besself(...,'z')
+## return a discrete space (Z) filter, W must be less than 1.
+##
+## [a,b,c,d] = besself(...)
+## return state-space matrices
+##
+## References:
+##
+## Proakis & Manolakis (1992). Digital Signal Processing. New York:
+## Macmillan Publishing Company.
+
+## Author: Paul Kienzle <pki...@us...>
+## Modified by: Doug Stewart <da...@sy...> Feb, 2003
+
+function [a, b, c, d] = besself (n, W, varargin)
+
+ if (nargin>4 || nargin<2) || (nargout>4 || nargout<2)
+ usage ("[b, a] or [z, p, g] or [a,b,c,d] = besself (n, W [, 'ftype'][,'z'])");
+ end
+
+ ## interpret the input parameters
+ if (!(length(n)==1 && n == round(n) && n > 0))
+ error ("besself: filter order n must be a positive integer");
+ end
+
+ stop = 0;
+ digital = 0;
+ for i=1:length(varargin)
+ switch varargin{i}
+ case 's', digital = 0;
+ case 'z', digital = 1;
+ case { 'high', 'stop' }, stop = 1;
+ case { 'low', 'pass' }, stop = 0;
+ otherwise, error ("besself: expected [high|stop] or [s|z]");
+ endswitch
+ endfor
+
+
+ [r, c]=size(W);
+ if (!(length(W)<=2 && (r==1 || c==1)))
+ error ("besself: frequency must be given as w0 or [w0, w1]");
+ elseif (!(length(W)==1 || length(W) == 2))
+ error ("besself: only one filter band allowed");
+ elseif (length(W)==2 && !(W(1) < W(2)))
+ error ("besself: first band edge must be smaller than second");
+ endif
+
+ if ( digital && !all(W >= 0 & W <= 1))
+ error ("besself: critical frequencies must be in (0 1)");
+ elseif ( !digital && !all(W >= 0 ))
+ error ("besself: critical frequencies must be in (0 inf)");
+ endif
+
+ ## Prewarp to the band edges to s plane
+ if digital
+ T = 2; # sampling frequency of 2 Hz
+ W = 2/T*tan(pi*W/T);
+ endif
+
+ ## Generate splane poles for the prototype bessel filter
+ [zero, pole, gain] = besselap(n);
+
+ ## splane frequency transform
+ [zero, pole, gain] = sftrans(zero, pole, gain, W, stop);
+
+ ## Use bilinear transform to convert poles to the z plane
+ if digital
+ [zero, pole, gain] = bilinear(zero, pole, gain, T);
+ endif
+
+ ## convert to the correct output form
+ if nargout==2,
+ a = real(gain*poly(zero));
+ b = real(poly(pole));
+ elseif nargout==3,
+ a = zero;
+ b = pole;
+ c = gain;
+ else
+ ## output ss results
+ [a, b, c, d] = zp2ss (zero, pole, gain);
+ endif
+
+endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <ha...@us...> - 2010-11-27 15:46:41
|
Revision: 7957
http://octave.svn.sourceforge.net/octave/?rev=7957&view=rev
Author: hauberg
Date: 2010-11-27 15:46:34 +0000 (Sat, 27 Nov 2010)
Log Message:
-----------
Several updates to 'invfreq' (from Pascal Depuis)
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/invfreq.m
trunk/octave-forge/main/signal/inst/invfreqs.m
trunk/octave-forge/main/signal/inst/invfreqz.m
Modified: trunk/octave-forge/main/signal/inst/invfreq.m
===================================================================
--- trunk/octave-forge/main/signal/inst/invfreq.m 2010-11-26 22:07:52 UTC (rev 7956)
+++ trunk/octave-forge/main/signal/inst/invfreq.m 2010-11-27 15:46:34 UTC (rev 7957)
@@ -61,27 +61,63 @@
% *final debugging
% 2007-08-03 Rolf Schirmacher <Rol...@Mu...>
% *replaced == by strcmp() for character string comparison
+% 2010-10-04 Pascal Dupuis <Pas...@uc...>
+% *added the ability to specify zeroes in B
% TODO: implement Steiglitz-McBride iterations
% TODO: improve numerical stability for high order filters (matlab is a bit better)
% TODO: modify to accept more argument configurations
-function [B,A] = invfreq(H,F,nB,nA,W,iter,tol,tr,plane)
- n = max(nA,nB);
+function [B, A, SigN] = invfreq(H, F, nB, nA, W, iter, tol, tr, plane, varargin)
+ if length(nB) > 1, zB = nB(2); nB = nB(1); else zB = 0; end
+ n = max(nA, nB);
m = n+1; mA = nA+1; mB = nB+1;
nF = length(F);
if nF ~= length(H), disp('invfreqz: length of H and F must be the same'); end;
- if nargin < 5, W = ones(1,nF); end;
+ if nargin < 5 || isempty(W), W = ones(1, nF); end;
if nargin < 6, iter = []; end
if nargin < 7 tol = []; end
- if nargin < 8, tr = ''; end
+ if nargin < 8 || isempty(tr), tr = ''; end
if nargin < 9, plane = 'z'; end
+ if nargin < 10, varargin = {}; end
if iter~=[], disp('no implementation for iter yet'),end
if tol ~=[], disp('no implementation for tol yet'),end
- if plane ~= 'z' & plane ~= 's', disp('invfreqz: Error in plane argument'), end
+ if (plane ~= 'z' && plane ~= 's'), disp('invfreqz: Error in plane argument'), end
- Ruu = zeros(mB,mB); Ryy = zeros(nA,nA); Ryu = zeros(nA,mB);
- Pu = zeros(mB,1); Py = zeros(nA,1);
+ [reg, prop ] = parseparams(varargin);
+ %# should we normalise freqs to avoid matrices with rank deficiency ?
+ norm = false;
+ %# by default, use Ordinary Least Square to solve normal equations
+ method = 'LS';
+ if length(prop) > 0
+ indi = 1; while indi <= length(prop)
+ switch prop{indi}
+ case 'norm'
+ if indi < length(prop) && ~ischar(prop{indi+1}),
+ norm = logical(prop{indi+1});
+ prop(indi:indi+1) = [];
+ continue
+ else
+ norm = true; prop(indi) = [];
+ continue
+ end
+ case 'method'
+ if indi < length(prop) && ischar(prop{indi+1}),
+ method = prop{indi+1};
+ prop(indi:indi+1) = [];
+ continue
+ else
+ error('invfreq.m: incorrect/missing method argument');
+ end
+ otherwise %# FIXME: just skip it for now
+ disp(sprintf("Ignoring unkown argument %s", varargin{indi}));
+ indi = indi + 1;
+ end
+ end
+ end
+
+ Ruu = zeros(mB, mB); Ryy = zeros(nA, nA); Ryu = zeros(nA, mB);
+ Pu = zeros(mB, 1); Py = zeros(nA,1);
if strcmp(tr,'trace')
disp(' ')
disp('Computing nonuniformly sampled, equation-error, rational filter.');
@@ -90,51 +126,121 @@
end
s = sqrt(-1)*F;
- if strcmp(plane,'z')
- if max(F)>pi || min(F)<0
+ switch plane
+ case 'z'
+ if max(F) > pi || min(F) < 0
disp('hey, you frequency is outside the range 0 to pi, making my own')
- F = linspace(0,pi,length(H));
+ F = linspace(0, pi, length(H));
s = sqrt(-1)*F;
end
s = exp(-s);
- end;
+ case 's'
+ if max(F) > 1e6 && n > 5,
+ if ~norm,
+ disp('Be carefull, there are risks of generating singular matrices');
+ disp('Call invfreqs as (..., "norm", true) to avoid it');
+ else
+ Fmax = max(F); s = sqrt(-1)*F/Fmax;
+ end
+ end
+ end
- for k=1:nF
+ for k=1:nF,
Zk = (s(k).^[0:n]).';
Hk = H(k);
aHks = Hk*conj(Hk);
Rk = (W(k)*Zk)*Zk';
rRk = real(Rk);
- Ruu = Ruu+rRk(1:mB,1:mB);
- Ryy = Ryy+aHks*rRk(2:mA,2:mA);
- Ryu = Ryu+real(Hk*Rk(2:mA,1:mB));
- Pu = Pu+W(k)*real(conj(Hk)*Zk(1:mB));
- Py = Py+(W(k)*aHks)*real(Zk(2:mA));
+ Ruu = Ruu + rRk(1:mB, 1:mB);
+ Ryy = Ryy + aHks*rRk(2:mA, 2:mA);
+ Ryu = Ryu + real(Hk*Rk(2:mA, 1:mB));
+ Pu = Pu + W(k)*real(conj(Hk)*Zk(1:mB));
+ Py = Py + (W(k)*aHks)*real(Zk(2:mA));
end;
+ Rr = ones(length(s), mB+nA); Zk = s;
+ for k = 1:min(nA, nB),
+ Rr(:, 1+k) = Zk;
+ Rr(:, mB+k) = -Zk.*H;
+ Zk = Zk.*s;
+ end
+ for k = 1+min(nA, nB):max(nA, nB)-1,
+ if k <= nB, Rr(:, 1+k) = Zk; end
+ if k <= nA, Rr(:, mB+k) = -Zk.*H; end
+ Zk = Zk.*s;
+ end
+ k = k+1;
+ if k <= nB, Rr(:, 1+k) = Zk; end
+ if k <= nA, Rr(:, mB+k) = -Zk.*H; end
+
+ %# complex to real equation system -- this ensures real solution
+ Rr = Rr(:, 1+zB:end);
+ Rr = [real(Rr); imag(Rr)]; Pr = [real(H(:)); imag(H(:))];
+ %# normal equations -- keep for ref
+ %# Rn= [Ruu(1+zB:mB, 1+zB:mB), -Ryu(:, 1+zB:mB)'; -Ryu(:, 1+zB:mB), Ryy];
+ %# Pn= [Pu(1+zB:mB); -Py];
- R = [Ruu,-Ryu';-Ryu,Ryy];
- P = [Pu;-Py];
- Theta = R\P;
+ switch method
+ case {'ls' 'LS'}
+ %# avoid scaling errors with Theta = R\P;
+ %# [Q, R] = qr([Rn Pn]); Theta = R(1:end, 1:end-1)\R(1:end, end);
+ [Q, R] = qr([Rr Pr], 0); Theta = R(1:end-1, 1:end-1)\R(1:end-1, end);
+ %# SigN = R(end, end-1);
+ SigN = R(end, end);
+ case {'tls' 'TLS'}
+ % [U, S, V] = svd([Rn Pn]);
+ % SigN = S(end, end-1);
+ % Theta = -V(1:end-1, end)/V(end, end);
+ [U, S, V] = svd([Rr Pr], 0);
+ SigN = S(end, end);
+ Theta = -V(1:end-1, end)/V(end, end);
+ case {'mls' 'MLS' 'qr' 'QR'}
+ % [Q, R] = qr([Rn Pn], 0);
+ %# solve the noised part -- DO NOT USE ECONOMY SIZE !
+ % [U, S, V] = svd(R(nA+1:end, nA+1:end));
+ % SigN = S(end, end-1);
+ % Theta = -V(1:end-1, end)/V(end, end);
+ %# unnoised part -- remove B contribution and back-substitute
+ % Theta = [R(1:nA, 1:nA)\(R(1:nA, end) - R(1:nA, nA+1:end-1)*Theta)
+ % Theta];
+ %# solve the noised part -- economy size OK as #rows > #columns
+ [Q, R] = qr([Rr Pr], 0);
+ eB = mB-zB; sA = eB+1;
+ [U, S, V] = svd(R(sA:end, sA:end));
+ %# noised (A) coefficients
+ Theta = -V(1:end-1, end)/V(end, end);
+ %# unnoised (B) part -- remove A contribution and back-substitute
+ Theta = [R(1:eB, 1:eB)\(R(1:eB, end) - R(1:eB, sA:end-1)*Theta)
+ Theta];
+ SigN = S(end, end);
+ otherwise
+ error("invfreq: unknown method %s", method);
+ end
- B = Theta(1:mB)';
- A = [1 Theta(mB+1:mB+nA)'];
+ B = [zeros(zB, 1); Theta(1:mB-zB)].';
+ A = [1; Theta(mB-zB+(1:nA))].';
+
if strcmp(plane,'s')
B = B(mB:-1:1);
A = A(mA:-1:1);
+ if norm, %# Frequencies were normalised -- unscale coefficients
+ Zk = Fmax.^[n:-1:0].';
+ for k = nB:-1:1+zB, B(k) = B(k)/Zk(k); end
+ for k = nA:-1:1, A(k) = A(k)/Zk(k); end
+ end
end
endfunction
%!demo
-%! order = 12; % order of test filter
+%! order = 6; % order of test filter
%! fc = 1/2; % sampling rate / 4
%! n = 128; % frequency grid size
-%! [B,A] = butter(order,fc);
-%! [H,w] = freqz(B,A,n);
-%! [Bh,Ah] = invfreq(H,w,order,order);
-%! [Hh,wh] = freqz(Bh,Ah,n);
+%! [B, A] = butter(order,fc);
+%! [H, w] = freqz(B,A,n);
+%! [Bh, Ah] = invfreq(H,w,order,order);
+%! [Hh, wh] = freqz(Bh,Ah,n);
%! xlabel("Frequency (rad/sample)");
%! ylabel("Magnitude");
-%! plot(w,[abs(H);abs(Hh)])
+%! plot(w,[abs(H), abs(Hh)])
%! legend('Original','Measured');
%! err = norm(H-Hh);
%! disp(sprintf('L2 norm of frequency response error = %f',err));
Modified: trunk/octave-forge/main/signal/inst/invfreqs.m
===================================================================
--- trunk/octave-forge/main/signal/inst/invfreqs.m 2010-11-26 22:07:52 UTC (rev 7956)
+++ trunk/octave-forge/main/signal/inst/invfreqs.m 2010-11-27 15:46:34 UTC (rev 7957)
@@ -49,9 +49,11 @@
% TODO: check invfreq.m for todo's
-function [B,A] = invfreqs(H,F,nB,nA,W,iter,tol,tr)
+function [B, A, SigN] = invfreqs(H,F,nB,nA,W,iter,tol,tr, varargin)
-if nargin < 8
+if nargin < 9
+ varargin = {};
+ if nargin < 8
tr = '';
if nargin < 7
tol = [];
@@ -62,22 +64,35 @@
end
end
end
+ end
end
% now for the real work
-[B,A] = invfreq(H,F,nB,nA,W,iter,tol,tr,'s');
+[B, A, SigN] = invfreq(H, F,nB, nA, W, iter, tol, tr, 's', varargin{:});
endfunction
%!demo
%! B = [1/2 1];
+%! B = [1 0 0];
%! A = [1 1];
-%! w = linspace(0,4,128);
-%! H = freqs(B,A,w);
-%! [Bh,Ah] = invfreqs(H,w,1,1);
+%! %#A = [1 36 630 6930 51975 270270 945945 2027025 2027025]/2027025;
+%! A = [1 21 210 1260 4725 10395 10395]/10395;
+%! A = [1 6 15 15]/15;
+%! w = linspace(0, 8, 128);
+%! H0 = freqs(B, A, w);
+%! Nn = (randn(size(w))+j*randn(size(w)))/sqrt(2);
+%! order = length(A) - 1;
+%! [Bh, Ah, Sig0] = invfreqs(H0, w, [length(B)-1 2], length(A)-1);
%! Hh = freqs(Bh,Ah,w);
+%! [BLS, ALS, SigLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "LS");
+%! HLS = freqs(BLS, ALS, w);
+%! [BTLS, ATLS, SigTLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "TLS");
+%! HTLS = freqs(BTLS, ATLS, w);
+%! [BMLS, AMLS, SigMLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "QR");
+%! HMLS = freqs(BMLS, AMLS, w);
%! xlabel("Frequency (rad/sec)");
%! ylabel("Magnitude");
-%! plot(w,[abs(H);abs(Hh)])
+%! plot(w,[abs(H0); abs(Hh)])
%! legend('Original','Measured');
-%! err = norm(H-Hh);
+%! err = norm(H0-Hh);
%! disp(sprintf('L2 norm of frequency response error = %f',err));
Modified: trunk/octave-forge/main/signal/inst/invfreqz.m
===================================================================
--- trunk/octave-forge/main/signal/inst/invfreqz.m 2010-11-26 22:07:52 UTC (rev 7956)
+++ trunk/octave-forge/main/signal/inst/invfreqz.m 2010-11-27 15:46:34 UTC (rev 7957)
@@ -44,9 +44,11 @@
% TODO: check invfreq.m for todo's
-function [B,A] = invfreqz(H,F,nB,nA,W,iter,tol,tr)
+function [B, A, SigN] = invfreqz(H, F, nB, nA, W, iter, tol, tr, varargin)
-if nargin < 8
+if nargin < 9
+ varargin = {};
+ if nargin < 8
tr = '';
if nargin < 7
tol = [];
@@ -57,26 +59,35 @@
end
end
end
+ end
end
+% now for the real work
+[B, A, SigN] = invfreq(H, F, nB, nA, W, iter, tol, tr, 'z', varargin{:});
-% now for the real work
-[B,A] = invfreq(H,F,nB,nA,W,iter,tol,tr,'z');
endfunction
%!demo
-%! order = 12; % order of test filter
+%! order = 9; % order of test filter
+%! % going to 10 or above leads to numerical instabilities and large errors
%! fc = 1/2; % sampling rate / 4
%! n = 128; % frequency grid size
-%! [B,A] = butter(order,fc);
-%! [H,w] = freqz(B,A,n);
-%! [Bh,Ah] = invfreqz(H,w,order,order);
-%! [Hh,wh] = freqz(Bh,Ah,n);
+%! [B0, A0] = butter(order, fc);
+%! [H0, w] = freqz(B0, A0, n);
+%! Nn = (randn(size(w))+j*randn(size(w)))/sqrt(2);
+%! [Bh, Ah, Sig0] = invfreqz(H0, w, order, order);
+%! [Hh, wh] = freqz(Bh, Ah, n);
+%! [BLS, ALS, SigLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "LS");
+%! HLS = freqz(BLS, ALS, n);
+%! [BTLS, ATLS, SigTLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "TLS");
+%! HTLS = freqz(BTLS, ATLS, n);
+%! [BMLS, AMLS, SigMLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "QR");
+%! HMLS = freqz(BMLS, AMLS, n);
%! xlabel("Frequency (rad/sample)");
%! ylabel("Magnitude");
-%! plot(w,[abs(H);abs(Hh)])
+%! plot(w,[abs(H0) abs(Hh)])
%! legend('Original','Measured');
-%! err = norm(H-Hh);
+%! err = norm(H0-Hh);
%! disp(sprintf('L2 norm of frequency response error = %f',err));
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <car...@us...> - 2011-09-29 17:29:28
|
Revision: 8635
http://octave.svn.sourceforge.net/octave/?rev=8635&view=rev
Author: carandraug
Date: 2011-09-29 17:29:22 +0000 (Thu, 29 Sep 2011)
Log Message:
-----------
ellip: moved ellipdemo.m to %!demo block inside ellip.m
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/ellip.m
Removed Paths:
-------------
trunk/octave-forge/main/signal/inst/ellipdemo.m
Modified: trunk/octave-forge/main/signal/inst/ellip.m
===================================================================
--- trunk/octave-forge/main/signal/inst/ellip.m 2011-09-29 17:24:45 UTC (rev 8634)
+++ trunk/octave-forge/main/signal/inst/ellip.m 2011-09-29 17:29:22 UTC (rev 8635)
@@ -133,3 +133,35 @@
endfunction
+%!demo
+%! clc
+%! disp('---------------------------> NELLIP 0.2 EXAMPLE <-------------------------')
+%! x=input("Let's calculate the filter order: [ENTER]");
+%! disp("")
+%! x=input("[n, Ws] = ellipord([.1 .2],.4,1,90); [ENTER]");
+%! [n, Ws] = ellipord([.1 .2],.4,1,90)
+%! disp("")
+%! x=input("Let's calculate the filter: [ENTER]");
+%! disp("")
+%! x=input("[b,a]=ellip(5,1,90,[.1,.2]); [ENTER]");
+%! [b,a]=ellip(5,1,90,[.1,.2])
+%! disp("")
+%! x=input("Let's calculate the frequency response: [ENTER]");
+%! disp("")
+%! x=input("[h,w]=freqz(b,a); [ENTER]");
+%! [h,w]=freqz(b,a);
+%!
+%! xlabel("Frequency");
+%! ylabel("abs(H[w])[dB]");
+%! axis([0,1,-100,0]);
+%! plot(w./pi, 20*log10(abs(h)), ';;')
+%!
+%! hold('on');
+%! x=ones(1,length(h));
+%! plot(w./pi, x.*-1, ';-1 dB;')
+%! plot(w./pi, x.*-90, ';-90 dB;')
+%! hold('off');
+%!
+%! xlabel("")
+%! ylabel("")
+%! clc
Deleted: trunk/octave-forge/main/signal/inst/ellipdemo.m
===================================================================
--- trunk/octave-forge/main/signal/inst/ellipdemo.m 2011-09-29 17:24:45 UTC (rev 8634)
+++ trunk/octave-forge/main/signal/inst/ellipdemo.m 2011-09-29 17:29:22 UTC (rev 8635)
@@ -1,46 +0,0 @@
-## Copyright (C) 2001 Paulo Neis <p_...@ya...>
-##
-## This program is free software; you can redistribute it and/or modify
-## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
-## (at your option) any later version.
-##
-## This program is distributed in the hope that it will be useful,
-## but WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-## GNU General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with this program; If not, see <http://www.gnu.org/licenses/>.
-
-clc
-disp('---------------------------> NELLIP 0.2 EXAMPLE <-------------------------')
-x=input("Let's calculate the filter order: [ENTER]");
-disp("")
-x=input("[n, Ws] = ellipord([.1 .2],.4,1,90); [ENTER]");
-[n, Ws] = ellipord([.1 .2],.4,1,90)
-disp("")
-x=input("Let's calculate the filter: [ENTER]");
-disp("")
-x=input("[b,a]=ellip(5,1,90,[.1,.2]); [ENTER]");
-[b,a]=ellip(5,1,90,[.1,.2])
-disp("")
-x=input("Let's calculate the frequency response: [ENTER]");
-disp("")
-x=input("[h,w]=freqz(b,a); [ENTER]");
-[h,w]=freqz(b,a);
-
-xlabel("Frequency");
-ylabel("abs(H[w])[dB]");
-axis([0,1,-100,0]);
-plot(w./pi, 20*log10(abs(h)), ';;')
-
-hold('on');
-x=ones(1,length(h));
-plot(w./pi, x.*-1, ';-1 dB;')
-plot(w./pi, x.*-90, ';-90 dB;')
-hold('off');
-
-xlabel("")
-ylabel("")
-clc
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <car...@us...> - 2011-09-30 13:42:51
|
Revision: 8640
http://octave.svn.sourceforge.net/octave/?rev=8640&view=rev
Author: carandraug
Date: 2011-09-30 13:42:44 +0000 (Fri, 30 Sep 2011)
Log Message:
-----------
impinvar/invimpinvar: new functions for signal package submitted by Rudy Eschauzier <res...@ya...>
Added Paths:
-----------
trunk/octave-forge/main/signal/inst/impinvar.m
trunk/octave-forge/main/signal/inst/invimpinvar.m
trunk/octave-forge/main/signal/inst/private/
trunk/octave-forge/main/signal/inst/private/h1_z_deriv.m
trunk/octave-forge/main/signal/inst/private/inv_residue.m
trunk/octave-forge/main/signal/inst/private/pad_poly.m
trunk/octave-forge/main/signal/inst/private/polyrev.m
trunk/octave-forge/main/signal/inst/private/to_real.m
Added: trunk/octave-forge/main/signal/inst/impinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/impinvar.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/impinvar.m 2011-09-30 13:42:44 UTC (rev 8640)
@@ -0,0 +1,103 @@
+## Copyright 2007 R.G.H. Eschauzier <res...@ya...>
+## Adapted by Carnë Draug on 2011 <car...@gm...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software
+## Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+## -*- texinfo -*-
+## @deftypefn{Function File} {[@var{b_out}, @var{a_out}] =} impinvar (@var{b}, @var{a}, @var{ts}, @var{tol})
+## @deftypefnx{Function File} {[@var{b_out}, @var{a_out}] =} impinvar (@var{b}, @var{a}, @var{ts})
+## @deftypefnx{Function File} {[@var{b_out}, @var{a_out}] =} impinvar (@var{b}, @var{a})
+## Converts analog filter with coefficients @var{b} and @var{a} to digital,
+## conserving impulse response.
+##
+## If @var{ts} is not specificied, or is an empty vector, it defaults to 1Hz.
+##
+## If @var{tol} is not specified, it defaults to 0.0001 (0.1%)
+## This function does the inverse of impinvar so that the following example should
+## restore the original values of @var{a} and @var{b}.
+##
+## @command{invimpinvar} implements the reverse of this function.
+## @example
+## [b, a] = impinvar (b, a);
+## [b, a] = invimpinvar (b, a);
+## @end example
+##
+## @seealso{bilinear, invimpinvar}
+## @end deftypefn
+
+function [b_out, a_out] = impinvar (b_in, a_in, ts = 1, tol = 0.0001)
+
+ if (nargin <2)
+ print_usage;
+ endif
+
+ ## to be compatible with the matlab implementation where an empty vector can
+ ## be used to get the default
+ if (isempty(ts))
+ ts = 1;
+ endif
+
+ [r_in, p_in, k_in] = residue(b_in, a_in); % partial fraction expansion
+
+ n = length(r_in); % Number of poles/residues
+
+ if (length(k_in)>0) % Greater than zero means we cannot do impulse invariance
+ error("Order numerator >= order denominator");
+ endif
+
+ r_out = zeros(1,n); % Residues of H(z)
+ p_out = zeros(1,n); % Poles of H(z)
+ k_out = 0; % Contstant term of H(z)
+
+ i=1;
+ while (i<=n)
+ m = 1;
+ first_pole = p_in(i); % Pole in the s-domain
+ while (i<n && abs(first_pole-p_in(i+1))<tol) % Multiple poles at p(i)
+ i++; % Next residue
+ m++; % Next multiplicity
+ endwhile
+ [r, p, k] = z_res(r_in(i-m+1:i), first_pole, ts); % Find z-domain residues
+ k_out += k; % Add direct term to output
+ p_out(i-m+1:i) = p; % Copy z-domain pole(s) to output
+ r_out(i-m+1:i) = r; % Copy z-domain residue(s) to output
+
+ i++; % Next s-domain residue/pole
+ endwhile
+
+ [b_out, a_out] = inv_residue(r_out, p_out, k_out);
+ a_out = to_real(a_out); % Get rid of spurious imaginary part
+ b_out = to_real(b_out);
+ b_out = polyreduce(b_out);
+
+endfunction
+
+## Convert residue vector for single and multiple poles in s-domain (located at sm) to
+## residue vector in z-domain. The variable k is the direct term of the result.
+function [r_out, p_out, k_out] = z_res (r_in, sm, ts);
+
+ p_out = exp(ts * sm); % z-domain pole
+ n = length(r_in); % Multiplicity of the pole
+ r_out = zeros(1,n); % Residue vector
+
+ %% First pole (no multiplicity)
+ k_out = r_in(1) * ts; % PFE of z/(z-p) = p/(z-p)+1; direct part
+ r_out(1) = r_in(1) * ts * p_out; % pole part of PFE
+
+ for i=(2:n) % Go through s-domain residues for multiple pole
+ r_out(1:i) += r_in(i) * polyrev(h1_z_deriv(i-1, p_out, ts)); % Add z-domain residues
+ endfor
+
+endfunction
Added: trunk/octave-forge/main/signal/inst/invimpinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/invimpinvar.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/invimpinvar.m 2011-09-30 13:42:44 UTC (rev 8640)
@@ -0,0 +1,98 @@
+## Copyright 2007 R.G.H. Eschauzier <res...@ya...>
+## Adapted by Carnë Draug on 2011 <car...@gm...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software
+## Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+## -*- texinfo -*-
+## @deftypefn{Function File} {[@var{b_out}, @var{a_out}] =} invimpinvar (@var{b}, @var{a}, @var{ts}, @var{tol})
+## @deftypefnx{Function File} {[@var{b_out}, @var{a_out}] =} invimpinvar (@var{b}, @var{a}, @var{ts})
+## @deftypefnx{Function File} {[@var{b_out}, @var{a_out}] =} invimpinvar (@var{b}, @var{a})
+## Converts digital filter with coefficients @var{b} and @var{a} to analog,
+## conserving impulse response.
+##
+## This function does the inverse of impinvar so that the following example should
+## restore the original values of @var{a} and @var{b}.
+## @example
+## [b, a] = impinvar (b, a);
+## [b, a] = invimpinvar (b, a);
+## @end example
+##
+## If @var{ts} is not specificied, or is an empty vector, it defaults to 1Hz.
+##
+## If @var{tol} is not specified, it defaults to 0.0001 (0.1%)
+##
+## @seealso{bilinear, impinvar}
+## @end deftypefn
+
+## Impulse invariant conversion from s to z domain
+function [b_out, a_out] = invimpinvar (b_in, a_in, ts = 1, tol = 0.0001)
+
+ if (nargin <2)
+ print_usage;
+ endif
+
+ [r_in, p_in, k_in] = residue(b_in, a_in); % partial fraction expansion
+
+ n = length(r_in); % Number of poles/residues
+
+ if (length(k_in) > 1) % Greater than one means we cannot do impulse invariance
+ error("Order numerator > order denominator");
+ endif
+
+ r_out = zeros(1,n); % Residues of H(s)
+ sm_out = zeros(1,n); % Poles of H(s)
+
+ i=1;
+ while (i<=n)
+ m=1;
+ first_pole = p_in(i); % Pole in the z-domain
+ while (i<n && abs(first_pole-p_in(i+1))<tol) % Multiple poles at p(i)
+ i++; % Next residue
+ m++; % Next multiplicity
+ endwhile
+ [r, sm, k] = inv_z_res(r_in(i-m+1:i), first_pole, ts); % Find s-domain residues
+ k_in -= k; % Just to check, should end up zero for physical system
+ sm_out(i-m+1:i) = sm; % Copy s-domain pole(s) to output
+ r_out(i-m+1:i) = r; % Copy s-domain residue(s) to output
+
+ i++; % Next z-domain residue/pole
+ endwhile
+ [b_out, a_out] = inv_residue(r_out, sm_out ,0);
+ a_out = to_real(a_out); % Get rid of spurious imaginary part
+ b_out = to_real(b_out);
+ b_out = polyreduce(b_out);
+
+endfunction
+
+## Inverse function of z_res (see impinvar source)
+function [r_out sm_out k_out] = inv_z_res (r_in,p_in,ts)
+
+ n = length(r_in); % multiplicity of the pole
+ r_in = r_in'; % From column vector to row vector
+
+ i=n;
+ while (i>1) % Go through residues starting from highest order down
+ r_out(i) = r_in(i) / ((ts * p_in)^i); % Back to binomial coefficient for highest order (always 1)
+ r_in(1:i) -= r_out(i) * polyrev(h1_z_deriv(i-1,p_in,ts)); % Subtract highest order result, leaving r_in(i) zero
+ i--;
+ endwhile
+
+ %% Single pole (no multiplicity)
+ r_out(1) = r_in(1) / ((ts * p_in));
+ k_out = r_in(1) / p_in;
+
+ sm_out = log(p_in) / ts;
+
+endfunction
Added: trunk/octave-forge/main/signal/inst/private/h1_z_deriv.m
===================================================================
--- trunk/octave-forge/main/signal/inst/private/h1_z_deriv.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/private/h1_z_deriv.m 2011-09-30 13:42:44 UTC (rev 8640)
@@ -0,0 +1,53 @@
+## Copyright 2007 R.G.H. Eschauzier <res...@ya...>
+## Adapted by Carnë Draug on 2011 <car...@gm...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software
+## Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+## This function is necessary for impinvar and invimpinvar of the signal package
+
+## Find {-zd/dz}^n*H1(z). I.e., first multiply by -z, then diffentiate, then multiply by -z, etc.
+## The result is (ts^(n+1))*(b(1)*p/(z-p)^1 + b(2)*p^2/(z-p)^2 + b(n+1)*p^(n+1)/(z-p)^(n+1)).
+## Works for n>0.
+function b = h1_z_deriv(n, p, ts)
+
+ %% Build the vector d that holds coefficients for all the derivatives of H1(z)
+ %% The results reads d(n)*z^(1)*(d/dz)^(1)*H1(z) + d(n-1)*z^(2)*(d/dz)^(2)*H1(z) +...+ d(1)*z^(n)*(d/dz)^(n)*H1(z)
+ d = (-1)^n; % Vector with the derivatives of H1(z)
+ for i=(1:n-1)
+ d = conv(d,[1 0]); % Shift result right (multiply by -z)
+ d += pad_poly(polyderiv(d),i+1); % Add the derivative
+ endfor
+
+ %% Build output vector
+ b = zeros(1,n+1);
+ for i=(1:n)
+ b += d(i) * pad_poly(h1_deriv(n-i+1, ts), n+1);
+ endfor
+
+ b *= ts^(n+1)/factorial(n);
+ for i=(1:n+1)
+ b(n-i+2) *= p^i; % Multiply coefficients with p^i, where i is the index of the coeff.
+ endfor
+
+endfunction
+
+## Find (z^n)*(d/dz)^n*H1(z), where H1(z)=ts*z/(z-p), ts=sampling period,
+## p=exp(sm*ts) and sm is the s-domain pole with multiplicity n+1.
+## The result is (ts^(n+1))*(b(1)*p/(z-p)^1 + b(2)*p^2/(z-p)^2 + b(n+1)*p^(n+1)/(z-p)^(n+1)),
+## where b(i) is the binomial coefficient bincoeff(n,i) times n!. Works for n>0.
+function b = h1_deriv(n,ts)
+ b = factorial(n)*bincoeff(n,0:n); % Binomial coefficients: [1], [1 1], [1 2 1], [1 3 3 1], etc.
+ b *= (-1)^n;
+endfunction
Added: trunk/octave-forge/main/signal/inst/private/inv_residue.m
===================================================================
--- trunk/octave-forge/main/signal/inst/private/inv_residue.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/private/inv_residue.m 2011-09-30 13:42:44 UTC (rev 8640)
@@ -0,0 +1,62 @@
+## Copyright 2007 R.G.H. Eschauzier <res...@ya...>
+## Adapted by Carnë Draug on 2011 <car...@gm...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software
+## Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+## This function is necessary for impinvar and invimpinvar of the signal package
+
+## Inverse of Octave residue function
+function [b_out, a_out] = inv_residue(r_in, p_in, k_in)
+
+ tol=0.0001;
+
+ n = length(r_in); % Number of poles/residues
+
+ k = 0; % Capture contstant term
+ if (length(k_in)==1) % A single direct term (order N = order D)
+ k = k_in(1); % Capture constant term
+ elseif (length(k_in)>1) % Greater than one means non-physical system
+ error("Order numerator > order denominator");
+ endif
+
+ a_out = 1; % Calculate denominator
+ for i=(1:n)
+ a_out = conv(a_out,[1 -p_in(i)]);
+ endfor
+
+ b_out = zeros(1,n+1);
+ b_out += k*a_out; % Constant term: add k times denominator to numerator
+ i=1;
+ while (i<=n)
+ term = [1 -p_in(i)]; % Term to be factored out
+ p = r_in(i)*deconv(a_out,term); % Residue times resulting polynomial
+ p = pad_poly(p,n+1); % Pad for proper length
+ b_out += p;
+
+ m = 1;
+ mterm = term;
+ first_pole = p_in(i);
+ while (i<n && abs(first_pole-p_in(i+1))<tol) % Multiple poles at p(i)
+ i++; %Next residue
+ m++;
+ mterm = conv(mterm,term); % Next multiplicity to be factored out
+ p = r_in(i)*deconv(a_out,mterm); % Resulting polynomial
+ p = pad_poly(p,n+1); % Pad for proper length
+ b_out += p;
+ endwhile
+ i++;
+ endwhile
+ b_out = polyreduce(b_out);
+endfunction
Added: trunk/octave-forge/main/signal/inst/private/pad_poly.m
===================================================================
--- trunk/octave-forge/main/signal/inst/private/pad_poly.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/private/pad_poly.m 2011-09-30 13:42:44 UTC (rev 8640)
@@ -0,0 +1,25 @@
+## Copyright 2007 R.G.H. Eschauzier <res...@ya...>
+## Adapted by Carnë Draug on 2011 <car...@gm...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software
+## Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+## This function is necessary for impinvar and invimpinvar of the signal package
+
+## Pad polynomial to length n
+function p_out = pad_poly(p_in,n)
+ l = length(p_in); % Length of input polynomial
+ p_out(n-l+1:n) = p_in(1:l); % Right shift for proper length
+ p_out(1:n-l) = 0; % Set first elements to zero
+endfunction
Added: trunk/octave-forge/main/signal/inst/private/polyrev.m
===================================================================
--- trunk/octave-forge/main/signal/inst/private/polyrev.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/private/polyrev.m 2011-09-30 13:42:44 UTC (rev 8640)
@@ -0,0 +1,26 @@
+## Copyright 2007 R.G.H. Eschauzier <res...@ya...>
+## Adapted by Carnë Draug on 2011 <car...@gm...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software
+## Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+## This function is necessary for impinvar and invimpinvar of the signal package
+
+## Reverse the coefficients of a polynomial
+function p_out = polyrev (p_in)
+ n = length(p_in);
+ for i=(1:n)
+ p_out(i) = p_in(n-i+1);
+ endfor
+endfunction
Added: trunk/octave-forge/main/signal/inst/private/to_real.m
===================================================================
--- trunk/octave-forge/main/signal/inst/private/to_real.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/private/to_real.m 2011-09-30 13:42:44 UTC (rev 8640)
@@ -0,0 +1,23 @@
+## Copyright 2007 R.G.H. Eschauzier <res...@ya...>
+## Adapted by Carnë Draug on 2011 <car...@gm...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software
+## Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+## This function is necessary for impinvar and invimpinvar of the signal package
+
+## Round complex number to nearest real number
+function p_out = to_real(p_in)
+ p_out = abs(p_in) .* sign(real(p_in));
+endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <car...@us...> - 2011-09-30 14:08:08
|
Revision: 8641
http://octave.svn.sourceforge.net/octave/?rev=8641&view=rev
Author: carandraug
Date: 2011-09-30 14:08:01 +0000 (Fri, 30 Sep 2011)
Log Message:
-----------
impinvar/invimpinvar: small bug fix where tol was hardcoded in auxiliary function
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/impinvar.m
trunk/octave-forge/main/signal/inst/invimpinvar.m
trunk/octave-forge/main/signal/inst/private/inv_residue.m
Modified: trunk/octave-forge/main/signal/inst/impinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/impinvar.m 2011-09-30 13:42:44 UTC (rev 8640)
+++ trunk/octave-forge/main/signal/inst/impinvar.m 2011-09-30 14:08:01 UTC (rev 8641)
@@ -77,7 +77,7 @@
i++; % Next s-domain residue/pole
endwhile
- [b_out, a_out] = inv_residue(r_out, p_out, k_out);
+ [b_out, a_out] = inv_residue(r_out, p_out, k_out, tol);
a_out = to_real(a_out); % Get rid of spurious imaginary part
b_out = to_real(b_out);
b_out = polyreduce(b_out);
Modified: trunk/octave-forge/main/signal/inst/invimpinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/invimpinvar.m 2011-09-30 13:42:44 UTC (rev 8640)
+++ trunk/octave-forge/main/signal/inst/invimpinvar.m 2011-09-30 14:08:01 UTC (rev 8641)
@@ -43,6 +43,12 @@
print_usage;
endif
+ ## to be compatible with the matlab implementation where an empty vector can
+ ## be used to get the default
+ if (isempty(ts))
+ ts = 1;
+ endif
+
[r_in, p_in, k_in] = residue(b_in, a_in); % partial fraction expansion
n = length(r_in); % Number of poles/residues
@@ -69,7 +75,7 @@
i++; % Next z-domain residue/pole
endwhile
- [b_out, a_out] = inv_residue(r_out, sm_out ,0);
+ [b_out, a_out] = inv_residue(r_out, sm_out , 0, tol);
a_out = to_real(a_out); % Get rid of spurious imaginary part
b_out = to_real(b_out);
b_out = polyreduce(b_out);
Modified: trunk/octave-forge/main/signal/inst/private/inv_residue.m
===================================================================
--- trunk/octave-forge/main/signal/inst/private/inv_residue.m 2011-09-30 13:42:44 UTC (rev 8640)
+++ trunk/octave-forge/main/signal/inst/private/inv_residue.m 2011-09-30 14:08:01 UTC (rev 8641)
@@ -18,10 +18,8 @@
## This function is necessary for impinvar and invimpinvar of the signal package
## Inverse of Octave residue function
-function [b_out, a_out] = inv_residue(r_in, p_in, k_in)
+function [b_out, a_out] = inv_residue(r_in, p_in, k_in, tol)
- tol=0.0001;
-
n = length(r_in); % Number of poles/residues
k = 0; % Capture contstant term
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <car...@us...> - 2011-09-30 14:17:58
|
Revision: 8642
http://octave.svn.sourceforge.net/octave/?rev=8642&view=rev
Author: carandraug
Date: 2011-09-30 14:17:52 +0000 (Fri, 30 Sep 2011)
Log Message:
-----------
impinvar/invimpinvar: change API to use sampling frequency instead of sampling time to be compatible with MatLab
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/impinvar.m
trunk/octave-forge/main/signal/inst/invimpinvar.m
Modified: trunk/octave-forge/main/signal/inst/impinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/impinvar.m 2011-09-30 14:08:01 UTC (rev 8641)
+++ trunk/octave-forge/main/signal/inst/impinvar.m 2011-09-30 14:17:52 UTC (rev 8642)
@@ -16,13 +16,13 @@
## Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
## -*- texinfo -*-
-## @deftypefn{Function File} {[@var{b_out}, @var{a_out}] =} impinvar (@var{b}, @var{a}, @var{ts}, @var{tol})
-## @deftypefnx{Function File} {[@var{b_out}, @var{a_out}] =} impinvar (@var{b}, @var{a}, @var{ts})
+## @deftypefn{Function File} {[@var{b_out}, @var{a_out}] =} impinvar (@var{b}, @var{a}, @var{fs}, @var{tol})
+## @deftypefnx{Function File} {[@var{b_out}, @var{a_out}] =} impinvar (@var{b}, @var{a}, @var{fs})
## @deftypefnx{Function File} {[@var{b_out}, @var{a_out}] =} impinvar (@var{b}, @var{a})
## Converts analog filter with coefficients @var{b} and @var{a} to digital,
## conserving impulse response.
##
-## If @var{ts} is not specificied, or is an empty vector, it defaults to 1Hz.
+## If @var{fs} is not specificied, or is an empty vector, it defaults to 1Hz.
##
## If @var{tol} is not specified, it defaults to 0.0001 (0.1%)
## This function does the inverse of impinvar so that the following example should
@@ -47,6 +47,8 @@
## be used to get the default
if (isempty(ts))
ts = 1;
+ else
+ ts = 1/ts; # we should be using sampling frequencies to be compatible with Matlab
endif
[r_in, p_in, k_in] = residue(b_in, a_in); % partial fraction expansion
Modified: trunk/octave-forge/main/signal/inst/invimpinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/invimpinvar.m 2011-09-30 14:08:01 UTC (rev 8641)
+++ trunk/octave-forge/main/signal/inst/invimpinvar.m 2011-09-30 14:17:52 UTC (rev 8642)
@@ -16,8 +16,8 @@
## Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
## -*- texinfo -*-
-## @deftypefn{Function File} {[@var{b_out}, @var{a_out}] =} invimpinvar (@var{b}, @var{a}, @var{ts}, @var{tol})
-## @deftypefnx{Function File} {[@var{b_out}, @var{a_out}] =} invimpinvar (@var{b}, @var{a}, @var{ts})
+## @deftypefn{Function File} {[@var{b_out}, @var{a_out}] =} invimpinvar (@var{b}, @var{a}, @var{fs}, @var{tol})
+## @deftypefnx{Function File} {[@var{b_out}, @var{a_out}] =} invimpinvar (@var{b}, @var{a}, @var{fs})
## @deftypefnx{Function File} {[@var{b_out}, @var{a_out}] =} invimpinvar (@var{b}, @var{a})
## Converts digital filter with coefficients @var{b} and @var{a} to analog,
## conserving impulse response.
@@ -29,7 +29,7 @@
## [b, a] = invimpinvar (b, a);
## @end example
##
-## If @var{ts} is not specificied, or is an empty vector, it defaults to 1Hz.
+## If @var{fs} is not specificied, or is an empty vector, it defaults to 1Hz.
##
## If @var{tol} is not specified, it defaults to 0.0001 (0.1%)
##
@@ -47,6 +47,8 @@
## be used to get the default
if (isempty(ts))
ts = 1;
+ else
+ ts = 1/ts; # we should be using sampling frequencies to be compatible with Matlab
endif
[r_in, p_in, k_in] = residue(b_in, a_in); % partial fraction expansion
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <car...@us...> - 2011-09-30 15:00:44
|
Revision: 8645
http://octave.svn.sourceforge.net/octave/?rev=8645&view=rev
Author: carandraug
Date: 2011-09-30 15:00:38 +0000 (Fri, 30 Sep 2011)
Log Message:
-----------
impinvar: bug fix by Juan Carbajal, of the bug fix to truncate values below tolerance
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/impinvar.m
trunk/octave-forge/main/signal/inst/invimpinvar.m
Modified: trunk/octave-forge/main/signal/inst/impinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/impinvar.m 2011-09-30 14:47:28 UTC (rev 8644)
+++ trunk/octave-forge/main/signal/inst/impinvar.m 2011-09-30 15:00:38 UTC (rev 8645)
@@ -1,5 +1,6 @@
## Copyright 2007 R.G.H. Eschauzier <res...@ya...>
-## Adapted by Carnë Draug on 2011 <car...@gm...>
+## Copyright (c) 2011 Carnë Draug <car...@gm...>
+## Copyright (c) 2011 Juan Pablo Carbajal <car...@if...>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
@@ -85,8 +86,8 @@
b_out = polyreduce(b_out);
## respect the required tolerance values
- a_out = a_out(find ((a_out < tol) == (a_out > -tol)));
- b_out = b_out(find ((b_out < tol) == (b_out > -tol)));
+ b_out(abs(b_out)<tol) = [];
+ a_out(abs(a_out)<tol) = [];
endfunction
Modified: trunk/octave-forge/main/signal/inst/invimpinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/invimpinvar.m 2011-09-30 14:47:28 UTC (rev 8644)
+++ trunk/octave-forge/main/signal/inst/invimpinvar.m 2011-09-30 15:00:38 UTC (rev 8645)
@@ -1,5 +1,5 @@
## Copyright 2007 R.G.H. Eschauzier <res...@ya...>
-## Adapted by Carnë Draug on 2011 <car...@gm...>
+## Copyright (c) 2011 Carnë Draug <car...@gm...>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <car...@us...> - 2011-10-01 20:33:02
|
Revision: 8648
http://octave.svn.sourceforge.net/octave/?rev=8648&view=rev
Author: carandraug
Date: 2011-10-01 20:32:55 +0000 (Sat, 01 Oct 2011)
Log Message:
-----------
impinvar/invimpinvar: removed private function pad_poly and replaced with octave-core prepad
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/impinvar.m
trunk/octave-forge/main/signal/inst/private/h1_z_deriv.m
trunk/octave-forge/main/signal/inst/private/inv_residue.m
Removed Paths:
-------------
trunk/octave-forge/main/signal/inst/private/pad_poly.m
Modified: trunk/octave-forge/main/signal/inst/impinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/impinvar.m 2011-09-30 19:06:11 UTC (rev 8647)
+++ trunk/octave-forge/main/signal/inst/impinvar.m 2011-10-01 20:32:55 UTC (rev 8648)
@@ -1,4 +1,4 @@
-## Copyright 2007 R.G.H. Eschauzier <res...@ya...>
+## Copyright (c) 2007 R.G.H. Eschauzier <res...@ya...>
## Copyright (c) 2011 Carnë Draug <car...@gm...>
## Copyright (c) 2011 Juan Pablo Carbajal <car...@if...>
##
Modified: trunk/octave-forge/main/signal/inst/private/h1_z_deriv.m
===================================================================
--- trunk/octave-forge/main/signal/inst/private/h1_z_deriv.m 2011-09-30 19:06:11 UTC (rev 8647)
+++ trunk/octave-forge/main/signal/inst/private/h1_z_deriv.m 2011-10-01 20:32:55 UTC (rev 8648)
@@ -27,13 +27,13 @@
d = (-1)^n; % Vector with the derivatives of H1(z)
for i=(1:n-1)
d = conv(d,[1 0]); % Shift result right (multiply by -z)
- d += pad_poly(polyderiv(d),i+1); % Add the derivative
+ d += prepad(polyderiv(d), i+1, 0) % Add the derivative
endfor
%% Build output vector
b = zeros(1,n+1);
for i=(1:n)
- b += d(i) * pad_poly(h1_deriv(n-i+1, ts), n+1);
+ b += d(i) * prepad(h1_deriv(n-i+1, ts), n+1, 0);
endfor
b *= ts^(n+1)/factorial(n);
Modified: trunk/octave-forge/main/signal/inst/private/inv_residue.m
===================================================================
--- trunk/octave-forge/main/signal/inst/private/inv_residue.m 2011-09-30 19:06:11 UTC (rev 8647)
+++ trunk/octave-forge/main/signal/inst/private/inv_residue.m 2011-10-01 20:32:55 UTC (rev 8648)
@@ -40,7 +40,7 @@
while (i<=n)
term = [1 -p_in(i)]; % Term to be factored out
p = r_in(i)*deconv(a_out,term); % Residue times resulting polynomial
- p = pad_poly(p,n+1); % Pad for proper length
+ p = prepad(p, n+1, 0); % Pad for proper length
b_out += p;
m = 1;
@@ -51,7 +51,7 @@
m++;
mterm = conv(mterm, term); % Next multiplicity to be factored out
p = r_in(i) * deconv(a_out, mterm); % Resulting polynomial
- p = pad_poly(p, n+1); % Pad for proper length
+ p = prepad(p, n+1, 0); % Pad for proper length
b_out += p;
endwhile
i++;
Deleted: trunk/octave-forge/main/signal/inst/private/pad_poly.m
===================================================================
--- trunk/octave-forge/main/signal/inst/private/pad_poly.m 2011-09-30 19:06:11 UTC (rev 8647)
+++ trunk/octave-forge/main/signal/inst/private/pad_poly.m 2011-10-01 20:32:55 UTC (rev 8648)
@@ -1,25 +0,0 @@
-## Copyright 2007 R.G.H. Eschauzier <res...@ya...>
-## Adapted by Carnë Draug on 2011 <car...@gm...>
-##
-## This program is free software; you can redistribute it and/or modify
-## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
-## (at your option) any later version.
-##
-## This program is distributed in the hope that it will be useful,
-## but WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-## GNU General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with this program; if not, write to the Free Software
-## Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-## This function is necessary for impinvar and invimpinvar of the signal package
-
-## Pad polynomial to length n
-function p_out = pad_poly(p_in,n)
- l = length(p_in); % Length of input polynomial
- p_out(n-l+1:n) = p_in(1:l); % Right shift for proper length
- p_out(1:n-l) = 0; % Set first elements to zero
-endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <car...@us...> - 2011-10-01 20:49:53
|
Revision: 8651
http://octave.svn.sourceforge.net/octave/?rev=8651&view=rev
Author: carandraug
Date: 2011-10-01 20:49:47 +0000 (Sat, 01 Oct 2011)
Log Message:
-----------
impinvar/invimpinvar: added reference to Thomas J. Cavicchi (1996)
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/impinvar.m
trunk/octave-forge/main/signal/inst/invimpinvar.m
Modified: trunk/octave-forge/main/signal/inst/impinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/impinvar.m 2011-10-01 20:45:27 UTC (rev 8650)
+++ trunk/octave-forge/main/signal/inst/impinvar.m 2011-10-01 20:49:47 UTC (rev 8651)
@@ -35,6 +35,9 @@
## [b, a] = invimpinvar (b, a);
## @end example
##
+## Reference: Thomas J. Cavicchi (1996) ``Impulse invariance and multiple-order
+## poles''. IEEE transactions on signal processing, Vol 40 (9): 2344--2347
+##
## @seealso{bilinear, invimpinvar}
## @end deftypefn
Modified: trunk/octave-forge/main/signal/inst/invimpinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/invimpinvar.m 2011-10-01 20:45:27 UTC (rev 8650)
+++ trunk/octave-forge/main/signal/inst/invimpinvar.m 2011-10-01 20:49:47 UTC (rev 8651)
@@ -33,6 +33,9 @@
##
## If @var{tol} is not specified, it defaults to 0.0001 (0.1%)
##
+## Reference: Thomas J. Cavicchi (1996) ``Impulse invariance and multiple-order
+## poles''. IEEE transactions on signal processing, Vol 40 (9): 2344--2347
+##
## @seealso{bilinear, impinvar}
## @end deftypefn
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <car...@us...> - 2011-10-14 17:13:28
|
Revision: 8746
http://octave.svn.sourceforge.net/octave/?rev=8746&view=rev
Author: carandraug
Date: 2011-10-14 17:13:21 +0000 (Fri, 14 Oct 2011)
Log Message:
-----------
impinvar: remove truncation of values below tolerance as it applies to the positions of the poles, but not necessarily to the coefficients of the z and s-domain polynomials. This becomes a problem for relatively high sample rates in impinvar() and low sample rates in invimpinvar().
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/impinvar.m
trunk/octave-forge/main/signal/inst/invimpinvar.m
Modified: trunk/octave-forge/main/signal/inst/impinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/impinvar.m 2011-10-13 19:18:05 UTC (rev 8745)
+++ trunk/octave-forge/main/signal/inst/impinvar.m 2011-10-14 17:13:21 UTC (rev 8746)
@@ -90,10 +90,6 @@
% Shift results right to account for calculating in z instead of z^-1
b_out(end)=[];
- ## respect the required tolerance values
- b_out(abs(b_out)<tol) = 0;
- a_out(abs(a_out)<tol) = 0;
-
endfunction
## Convert residue vector for single and multiple poles in s-domain (located at sm) to
@@ -118,7 +114,7 @@
%!function err = stozerr(bs,as,fs)
%!
%! % number of time steps
-%! n=10;
+%! n=100;
%!
%! % impulse invariant transform to z-domain
%! [bz az]=impinvar(bs,as,fs);
@@ -137,15 +133,15 @@
%! err=sqrt(sum((yz*fs.-ys).^2)/length(ys));
%! endfunction
%!
-%!assert(stozerr([1],[1 1],0.1),0,0.0001);
-%!assert(stozerr([1],[1 2 1],0.1),0,0.0001);
-%!assert(stozerr([1 1],[1 2 1],0.1),0,0.0001);
-%!assert(stozerr([1],[1 3 3 1],0.1),0,0.0001);
-%!assert(stozerr([1 1],[1 3 3 1],0.1),0,0.0001);
-%!assert(stozerr([1 1 1],[1 3 3 1],0.1),0,0.0001);
-%!assert(stozerr([1],[1 0 1],0.1),0,0.0001);
-%!assert(stozerr([1 1],[1 0 1],0.1),0,0.0001);
-%!assert(stozerr([1],[1 0 2 0 1],0.1),0,0.0001);
-%!assert(stozerr([1 1],[1 0 2 0 1],0.1),0,0.0001);
-%!assert(stozerr([1 1 1],[1 0 2 0 1],0.1),0,0.0001);
-%!assert(stozerr([1 1 1 1],[1 0 2 0 1],0.1),0,0.0001);
+%!assert(stozerr([1],[1 1],100),0,0.0001);
+%!assert(stozerr([1],[1 2 1],100),0,0.0001);
+%!assert(stozerr([1 1],[1 2 1],100),0,0.0001);
+%!assert(stozerr([1],[1 3 3 1],100),0,0.0001);
+%!assert(stozerr([1 1],[1 3 3 1],100),0,0.0001);
+%!assert(stozerr([1 1 1],[1 3 3 1],100),0,0.0001);
+%!assert(stozerr([1],[1 0 1],100),0,0.0001);
+%!assert(stozerr([1 1],[1 0 1],100),0,0.0001);
+%!assert(stozerr([1],[1 0 2 0 1],100),0,0.0001);
+%!assert(stozerr([1 1],[1 0 2 0 1],100),0,0.0001);
+%!assert(stozerr([1 1 1],[1 0 2 0 1],100),0,0.0001);
+%!assert(stozerr([1 1 1 1],[1 0 2 0 1],100),0,0.0001);
Modified: trunk/octave-forge/main/signal/inst/invimpinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/invimpinvar.m 2011-10-13 19:18:05 UTC (rev 8745)
+++ trunk/octave-forge/main/signal/inst/invimpinvar.m 2011-10-14 17:13:21 UTC (rev 8746)
@@ -1,4 +1,4 @@
-## Copyright 2007 R.G.H. Eschauzier <res...@ya...>
+## Copyright (c) 2007 R.G.H. Eschauzier <res...@ya...>
## Copyright (c) 2011 Carnë Draug <car...@gm...>
##
## This program is free software; you can redistribute it and/or modify
@@ -86,10 +86,6 @@
a_out = to_real(a_out); % Get rid of spurious imaginary part
b_out = to_real(b_out);
- ## respect the required tolerance values
- b_out(abs(b_out)<tol) = 0;
- a_out(abs(a_out)<tol) = 0;
-
b_out = polyreduce(b_out);
endfunction
@@ -118,7 +114,7 @@
%!function err = ztoserr(bz,az,fs)
%!
%! % number of time steps
-%! n=10;
+%! n=100;
%!
%! % make sure system is realizable (no delays)
%! bz=prepad(bz,length(az)-1,0,2);
@@ -140,15 +136,15 @@
%! err=sqrt(sum((yz*fs.-ys).^2)/length(ys));
%! endfunction
%!
-%!assert(ztoserr([1],[1 -0.5],10),0,0.0001);
-%!assert(ztoserr([1],[1 -1 0.25],10),0,0.0001);
-%!assert(ztoserr([1 1],[1 -1 0.25],10),0,0.0001);
-%!assert(ztoserr([1],[1 -1.5 0.75 -0.125],10),0,0.0001);
-%!assert(ztoserr([1 1],[1 -1.5 0.75 -0.125],10),0,0.0001);
-%!assert(ztoserr([1 1 1],[1 -1.5 0.75 -0.125],10),0,0.0001);
-%!assert(ztoserr([1],[1 0 0.25],10),0,0.0001);
-%!assert(ztoserr([1 1],[1 0 0.25],10),0,0.0001);
-%!assert(ztoserr([1],[1 0 0.5 0 0.0625],10),0,0.0001);
-%!assert(ztoserr([1 1],[1 0 0.5 0 0.0625],10),0,0.0001);
-%!assert(ztoserr([1 1 1],[1 0 0.5 0 0.0625],10),0,0.0001);
-%!assert(ztoserr([1 1 1 1],[1 0 0.5 0 0.0625],10),0,0.0001);
+%!assert(ztoserr([1],[1 -0.5],0.01),0,0.0001);
+%!assert(ztoserr([1],[1 -1 0.25],0.01),0,0.0001);
+%!assert(ztoserr([1 1],[1 -1 0.25],0.01),0,0.0001);
+%!assert(ztoserr([1],[1 -1.5 0.75 -0.125],0.01),0,0.0001);
+%!assert(ztoserr([1 1],[1 -1.5 0.75 -0.125],0.01),0,0.0001);
+%!assert(ztoserr([1 1 1],[1 -1.5 0.75 -0.125],0.01),0,0.0001);
+%!assert(ztoserr([1],[1 0 0.25],0.01),0,0.0001);
+%!assert(ztoserr([1 1],[1 0 0.25],0.01),0,0.0001);
+%!assert(ztoserr([1],[1 0 0.5 0 0.0625],0.01),0,0.0001);
+%!assert(ztoserr([1 1],[1 0 0.5 0 0.0625],0.01),0,0.0001);
+%!assert(ztoserr([1 1 1],[1 0 0.5 0 0.0625],0.01),0,0.0001);
+%!assert(ztoserr([1 1 1 1],[1 0 0.5 0 0.0625],0.01),0,0.0001);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2011-11-04 19:15:30
|
Revision: 8984
http://octave.svn.sourceforge.net/octave/?rev=8984&view=rev
Author: jpicarbajal
Date: 2011-11-04 19:15:24 +0000 (Fri, 04 Nov 2011)
Log Message:
-----------
signal. Fixing docstring and Matlab-style short-circuit operator warnings
Before release 1.1.0
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/data2fun.m
trunk/octave-forge/main/signal/inst/firls.m
trunk/octave-forge/main/signal/inst/fracshift.m
trunk/octave-forge/main/signal/inst/resample.m
trunk/octave-forge/main/signal/inst/sos2tf.m
Modified: trunk/octave-forge/main/signal/inst/data2fun.m
===================================================================
--- trunk/octave-forge/main/signal/inst/data2fun.m 2011-11-04 19:07:00 UTC (rev 8983)
+++ trunk/octave-forge/main/signal/inst/data2fun.m 2011-11-04 19:15:24 UTC (rev 8984)
@@ -38,6 +38,8 @@
%% @item interp
%% Type of interpolation. See @code{interp1}.
%%
+%% @end table
+%%
%% @seealso{interp1}
%% @end deftypefn
Modified: trunk/octave-forge/main/signal/inst/firls.m
===================================================================
--- trunk/octave-forge/main/signal/inst/firls.m 2011-11-04 19:07:00 UTC (rev 8983)
+++ trunk/octave-forge/main/signal/inst/firls.m 2011-11-04 19:15:24 UTC (rev 8984)
@@ -39,7 +39,7 @@
function coef = firls(N, frequencies, pass, weight, str);
- if nargin<3 | nargin>6
+ if nargin<3 || nargin>6
usage("");
end
if nargin==3
Modified: trunk/octave-forge/main/signal/inst/fracshift.m
===================================================================
--- trunk/octave-forge/main/signal/inst/fracshift.m 2011-11-04 19:07:00 UTC (rev 8983)
+++ trunk/octave-forge/main/signal/inst/fracshift.m 2011-11-04 19:15:24 UTC (rev 8984)
@@ -21,12 +21,12 @@
## @end deftypefn
## @seealso{circshift}
-## Ref [1] A. V. Oppenheim, R. W. Schafer and J. R. Buck,
+## Ref [1] A. V. Oppenheim, R. W. Schafer and J. R. Buck,
## Discrete-time signal processing, Signal processing series,
## Prentice-Hall, 1999
##
## Ref [2] T.I. Laakso, V. Valimaki, M. Karjalainen and U.K. Laine
-## Splitting the unit delay, IEEE Signal Processing Magazine,
+## Splitting the unit delay, IEEE Signal Processing Magazine,
## vol. 13, no. 1, pp 30--59 Jan 1996
function [y, h] = fracshift( x, d, h )
@@ -47,35 +47,35 @@
if (nargin < 4)
## properties of the interpolation filter
-
+
log10_rejection = -3.0;
stopband_cutoff_f = 1.0 / 2.0;
roll_off_width = stopband_cutoff_f / 10;
-
+
## determine filter length
## use empirical formula from [1] Chap 7, Eq. (7.63) p 476
-
+
rejection_dB = -20.0*log10_rejection;
L = ceil((rejection_dB-8.0) / (28.714 * roll_off_width));
-
+
## ideal sinc filter
-
+
t=(-L:L)';
- ideal_filter=2*stopband_cutoff_f*sinc(2*stopband_cutoff_f*(t-(d-fix(d))));
-
+ ideal_filter=2*stopband_cutoff_f*sinc(2*stopband_cutoff_f*(t-(d-fix(d))));
+
## determine parameter of Kaiser window
## use empirical formula from [1] Chap 7, Eq. (7.62) p 474
-
- if ((rejection_dB>=21)&(rejection_dB<=50))
+
+ if ((rejection_dB>=21) && (rejection_dB<=50))
beta = 0.5842 * (rejection_dB-21.0)^0.4 + 0.07886 * (rejection_dB-21.0);
elseif (rejection_dB>50)
beta = 0.1102 * (rejection_dB-8.7);
else
beta = 0.0;
endif
-
+
## apodize ideal (sincard) filter response
-
+
m = 2*L;
t = (0 : m)' - (d-fix(d));
t = 2 * beta / m * sqrt (t .* (m - t));
Modified: trunk/octave-forge/main/signal/inst/resample.m
===================================================================
--- trunk/octave-forge/main/signal/inst/resample.m 2011-11-04 19:07:00 UTC (rev 8983)
+++ trunk/octave-forge/main/signal/inst/resample.m 2011-11-04 19:15:24 UTC (rev 8984)
@@ -24,7 +24,7 @@
## Digital Signal Processing: Principles, Algorithms, and Applications,
## 4th ed., Prentice Hall, 2007. Chap. 6
##
-## Ref [2] A. V. Oppenheim, R. W. Schafer and J. R. Buck,
+## Ref [2] A. V. Oppenheim, R. W. Schafer and J. R. Buck,
## Discrete-time signal processing, Signal processing series,
## Prentice-Hall, 1999
## @end deftypefn
@@ -52,42 +52,42 @@
if (nargin < 4)
## properties of the antialiasing filter
-
+
log10_rejection = -3.0;
stopband_cutoff_f = 1.0/(2.0 * max(p,q));
roll_off_width = stopband_cutoff_f / 10.0;
-
+
## determine filter length
## use empirical formula from [2] Chap 7, Eq. (7.63) p 476
-
+
rejection_dB = -20.0*log10_rejection;
L = ceil((rejection_dB-8.0) / (28.714 * roll_off_width));
-
+
## ideal sinc filter
-
+
t=(-L:L)';
- ideal_filter=2*p*stopband_cutoff_f*sinc(2*stopband_cutoff_f*t);
-
+ ideal_filter=2*p*stopband_cutoff_f*sinc(2*stopband_cutoff_f*t);
+
## determine parameter of Kaiser window
## use empirical formula from [2] Chap 7, Eq. (7.62) p 474
-
- if ((rejection_dB>=21)&(rejection_dB<=50))
+
+ if ((rejection_dB>=21) && (rejection_dB<=50))
beta = 0.5842 * (rejection_dB-21.0)^0.4 + 0.07886 * (rejection_dB-21.0);
elseif (rejection_dB>50)
beta = 0.1102 * (rejection_dB-8.7);
else
beta = 0.0;
endif
-
+
## apodize ideal filter response
-
+
h=kaiser(2*L+1,beta).*ideal_filter;
-
+
endif
## check if input is a row vector
isrowvector=false;
-if ((rows(x)==1)&(columns(x)>1))
+if ((rows(x)==1) && (columns(x)>1))
x=x(:);
isrowvector=true;
endif
@@ -163,7 +163,7 @@
%! xx=sin(2*pi*f0/r*tt' + phi0);
%! t0=ceil((length(h)-1)/2/q);
%! idx=t0+1:NN-t0;
-%! reject(n+1)=max(abs(y(idx)));
+%! reject(n+1)=max(abs(y(idx)));
%! endfor;
%! rolloff=.1;
%! rejection=10^-3;
Modified: trunk/octave-forge/main/signal/inst/sos2tf.m
===================================================================
--- trunk/octave-forge/main/signal/inst/sos2tf.m 2011-11-04 19:07:00 UTC (rev 8983)
+++ trunk/octave-forge/main/signal/inst/sos2tf.m 2011-11-04 19:15:24 UTC (rev 8984)
@@ -23,7 +23,7 @@
%% @item
%% @var{sos} = matrix of series second-order sections, one per row:@*
%% @var{sos} = [@var{B1}.' @var{A1}.'; ...; @var{BN}.' @var{AN}.'], where@*
-%% @code{@var{B1}.'==[b0 b1 b2] and @var{A1}.'==[1 a1 a2]} for
+%% @code{@var{B1}.'==[b0 b1 b2] and @var{A1}.'==[1 a1 a2]} for
%% section 1, etc.@*
%% b0 must be nonzero for each section.@*
%% See @code{filter()} for documentation of the
@@ -33,31 +33,49 @@
%% @var{Bscale} is an overall gain factor that effectively scales
%% the output @var{B} vector (or any one of the input @var{B}i vectors).
%% @end itemize
-%%
+%%
%% RETURNED:
%% @var{B} and @var{A} are vectors specifying the digital filter @math{H(z) = B(z)/A(z)}.
%% See @code{filter()} for further details.
-%%
+%%
%% @seealso{tf2sos zp2sos sos2pz zp2tf tf2zp}
%% @end deftypefn
function [B,A] = sos2tf(sos,Bscale)
-if nargin<2, Bscale=1; end
-[N,M] = size(sos);
-if M~=6, error('sos2tf: sos matrix should be N by 6'); end
-A = 1;
-B = 1;
-for i=1:N
- B = conv(B, sos(i,1:3));
- A = conv(A, sos(i,4:6));
-end
-nB = length(B);
-while nB & B(nB)==0, B=B(1:nB-1); nB=length(B); end
-nA = length(A);
-while nA & A(nA)==0, A=A(1:nA-1); nA=length(A); end
-B = B * Bscale;
+ if nargin<2
+ Bscale=1;
+ end
+ [N,M] = size(sos);
+
+ if M~=6
+ error('sos2tf: sos matrix should be N by 6');
+ end
+
+ A = 1;
+ B = 1;
+
+ for i=1:N
+ B = conv(B, sos(i,1:3));
+ A = conv(A, sos(i,4:6));
+ end
+
+ nB = length(B);
+ while nB && B(nB)==0
+ B=B(1:nB-1);
+ nB=length(B);
+ end
+
+ nA = length(A);
+ while nA && A(nA)==0
+ A=A(1:nA-1);
+ nA=length(A);
+ end
+ B = B * Bscale;
+
+endfunction
+
%!test
%! B=[1 1];
%! A=[1 0.5];
@@ -66,7 +84,7 @@
%! assert({Bh,Ah},{B,A},10*eps);
%!test
-%! B=[1 0 0 0 0 1];
+%! B=[1 0 0 0 0 1];
%! A=[1 0 0 0 0 0.9];
%! [sos,g] = tf2sos(B,A);
%! [Bh,Ah] = sos2tf(sos,g);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <car...@us...> - 2011-11-05 13:45:57
|
Revision: 8994
http://octave.svn.sourceforge.net/octave/?rev=8994&view=rev
Author: carandraug
Date: 2011-11-05 13:45:50 +0000 (Sat, 05 Nov 2011)
Log Message:
-----------
signal package: added copyright word to public domain notice so that octave's get help don'i confuse it with help text block
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/downsample.m
trunk/octave-forge/main/signal/inst/dst.m
trunk/octave-forge/main/signal/inst/flattopwin.m
trunk/octave-forge/main/signal/inst/idst.m
trunk/octave-forge/main/signal/inst/upsample.m
Modified: trunk/octave-forge/main/signal/inst/downsample.m
===================================================================
--- trunk/octave-forge/main/signal/inst/downsample.m 2011-11-05 08:55:51 UTC (rev 8993)
+++ trunk/octave-forge/main/signal/inst/downsample.m 2011-11-05 13:45:50 UTC (rev 8994)
@@ -1,4 +1,4 @@
-## Author: Paul Kienzle <pki...@us...>
+## Copyright (C) 2007 Paul Kienzle <pki...@us...>
## This function is public domain
## -*- texinfo -*-
@@ -36,4 +36,3 @@
%!assert(downsample([1,2;3,4;5,6;7,8;9,10],2),[1,2;5,6;9,10]);
%!assert(downsample([1,2,3,4,5],2,1),[2,4]);
%!assert(downsample([1,2;3,4;5,6;7,8;9,10],2,1),[3,4;7,8]);
-
Modified: trunk/octave-forge/main/signal/inst/dst.m
===================================================================
--- trunk/octave-forge/main/signal/inst/dst.m 2011-11-05 08:55:51 UTC (rev 8993)
+++ trunk/octave-forge/main/signal/inst/dst.m 2011-11-05 13:45:50 UTC (rev 8994)
@@ -1,7 +1,5 @@
-## Author: Paul Kienzle <pki...@us...>
-## 2006-12-05
-## * initial release
-## This program is public domain
+## Copyright (C) 2006 Paul Kienzle <pki...@us...>
+## This function is public domain
## -*- texinfo -*-
## @deftypefn {Function File} @var{y} = dst (@var{x})
Modified: trunk/octave-forge/main/signal/inst/flattopwin.m
===================================================================
--- trunk/octave-forge/main/signal/inst/flattopwin.m 2011-11-05 08:55:51 UTC (rev 8993)
+++ trunk/octave-forge/main/signal/inst/flattopwin.m 2011-11-05 13:45:50 UTC (rev 8994)
@@ -1,4 +1,5 @@
-## This program is public domain.
+## Copyright (C) 2004 Paul Kienzle <pki...@us...>
+## This function is public domain
## flattopwin(n, [periodic|symmetric])
##
Modified: trunk/octave-forge/main/signal/inst/idst.m
===================================================================
--- trunk/octave-forge/main/signal/inst/idst.m 2011-11-05 08:55:51 UTC (rev 8993)
+++ trunk/octave-forge/main/signal/inst/idst.m 2011-11-05 13:45:50 UTC (rev 8994)
@@ -1,8 +1,5 @@
-## Author: Paul Kienzle <pki...@us...>
-## 2006-12-05
-## * initial release
-##
-## This program is public domain
+## Copyright (C) 2006 Paul Kienzle <pki...@us...>
+## This function is public domain
## -*- texinfo -*-
## @deftypefn {Function File} @var{y} = idst (@var{x})
@@ -13,6 +10,7 @@
## columns of the the matrix.
## @end deftypefn
## @seealso{dst}
+
function x = idst (y, n)
if (nargin < 1 || nargin > 2)
Modified: trunk/octave-forge/main/signal/inst/upsample.m
===================================================================
--- trunk/octave-forge/main/signal/inst/upsample.m 2011-11-05 08:55:51 UTC (rev 8993)
+++ trunk/octave-forge/main/signal/inst/upsample.m 2011-11-05 13:45:50 UTC (rev 8994)
@@ -1,4 +1,4 @@
-## Author: Paul Kienzle <pki...@us...>
+## Copyright (C) 2007 Paul Kienzle <pki...@us...>
## This function is public domain
## -*- texinfo -*-
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <car...@us...> - 2011-11-05 14:03:14
|
Revision: 8996
http://octave.svn.sourceforge.net/octave/?rev=8996&view=rev
Author: carandraug
Date: 2011-11-05 14:03:07 +0000 (Sat, 05 Nov 2011)
Log Message:
-----------
signal package:
* use print_usage
* use ... instead of \ as continuation line
* small fixes on help text
* set defaults in the function definition line
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/downsample.m
trunk/octave-forge/main/signal/inst/dst.m
trunk/octave-forge/main/signal/inst/idst.m
trunk/octave-forge/main/signal/inst/upsample.m
Modified: trunk/octave-forge/main/signal/inst/downsample.m
===================================================================
--- trunk/octave-forge/main/signal/inst/downsample.m 2011-11-05 14:02:06 UTC (rev 8995)
+++ trunk/octave-forge/main/signal/inst/downsample.m 2011-11-05 14:03:07 UTC (rev 8996)
@@ -2,26 +2,27 @@
## This function is public domain
## -*- texinfo -*-
-## @deftypefn {Function File} @var{y} = downsample(@var{x},@var{n})
+## @deftypefn {Function File} {@var{y} =} downsample (@var{x}, @var{n})
+## @deftypefnx {Function File} {@var{y} =} downsample (@var{x}, @var{n}, @var{offset})
## Downsample the signal, selecting every nth element. If @var{x}
## is a matrix, downsample every column.
##
-## For most signals you will want to use decimate() instead since
+## For most signals you will want to use @code{decimate} instead since
## it prefilters the high frequency components of the signal and
## avoids aliasing effects.
##
-## @deftypefnx {Function File} @var{y} = downsample(@var{x},@var{n},@var{offset})
-## Select every nth element starting at sample @var{offset}.
+## If @var{offset} is defined, select every nth element starting at
+## sample @var{offset}.
+## @seealso{decimate, interp, resample, upfirdn, upsample}
## @end deftypefn
-## @seealso{decimate, interp, resample, upfirdn, upsample}
-function y = downsample(x,n,phase)
- if nargin<2 || nargin>3, usage('downsample(x,n,[offset]'); end
- if nargin==2, phase = 0; end
+function y = downsample (x, n, phase = 0)
+ if nargin<2 || nargin>3, print_usage; end
+
if phase > n - 1
- warning("This is incompatible with Matlab (phase = 0:n-1). See\
- octave-forge signal package release notes for details." )
+ warning("This is incompatible with Matlab (phase = 0:n-1). See ...
+ octave-forge signal package release notes for details.")
end
if isvector(x)
Modified: trunk/octave-forge/main/signal/inst/dst.m
===================================================================
--- trunk/octave-forge/main/signal/inst/dst.m 2011-11-05 14:02:06 UTC (rev 8995)
+++ trunk/octave-forge/main/signal/inst/dst.m 2011-11-05 14:03:07 UTC (rev 8996)
@@ -2,8 +2,8 @@
## This function is public domain
## -*- texinfo -*-
-## @deftypefn {Function File} @var{y} = dst (@var{x})
-## @deftypefnx {Function File} @var{y} = dst (@var{x}, @var{n})
+## @deftypefn {Function File} {@var{y} =} dst (@var{x})
+## @deftypefnx {Function File} {@var{y} =} dst (@var{x}, @var{n})
## Computes the type I discrete sine transform of @var{x}. If @var{n} is given,
## then @var{x} is padded or trimmed to length @var{n} before computing the transform.
## If @var{x} is a matrix, compute the transform along the columns of the
@@ -11,17 +11,19 @@
##
## The discrete sine transform X of x can be defined as follows:
##
+## @verbatim
## N
## X[k] = sum x[n] sin (pi n k / (N+1) ), k = 1, ..., N
## n=1
+## @end verbatim
##
+## @seealso{idst}
## @end deftypefn
-## @seealso{idst}
function y = dst (x, n)
if (nargin < 1 || nargin > 2)
- usage ("y = dst(x [, n])");
+ print_usage;
endif
transpose = (rows (x) == 1);
@@ -48,10 +50,8 @@
if transpose, y = y.'; endif
-
endfunction
-
%!test
%! x = log(linspace(0.1,1,32));
%! y = dst(x);
Modified: trunk/octave-forge/main/signal/inst/idst.m
===================================================================
--- trunk/octave-forge/main/signal/inst/idst.m 2011-11-05 14:02:06 UTC (rev 8995)
+++ trunk/octave-forge/main/signal/inst/idst.m 2011-11-05 14:03:07 UTC (rev 8996)
@@ -2,19 +2,19 @@
## This function is public domain
## -*- texinfo -*-
-## @deftypefn {Function File} @var{y} = idst (@var{x})
-## @deftypefnx {Function File} @var{y} = idst (@var{x}, @var{n})
+## @deftypefn {Function File} {@var{y} =} idst (@var{x})
+## @deftypefnx {Function File} {@var{y} =} idst (@var{x}, @var{n})
## Computes the inverse type I discrete sine transform of @var{y}. If @var{n} is
## given, then @var{y} is padded or trimmed to length @var{n} before computing
## the transform. If @var{y} is a matrix, compute the transform along the
## columns of the the matrix.
+## @seealso{dst}
## @end deftypefn
-## @seealso{dst}
function x = idst (y, n)
if (nargin < 1 || nargin > 2)
- usage ("x = idst(y [, n])");
+ print_usage;
endif
if nargin == 1,
Modified: trunk/octave-forge/main/signal/inst/upsample.m
===================================================================
--- trunk/octave-forge/main/signal/inst/upsample.m 2011-11-05 14:02:06 UTC (rev 8995)
+++ trunk/octave-forge/main/signal/inst/upsample.m 2011-11-05 14:03:07 UTC (rev 8996)
@@ -2,21 +2,22 @@
## This function is public domain
## -*- texinfo -*-
-## @deftypefn {Function File} @var{y} = upsample(@var{x},@var{n})
-## Upsample the signal, inserting n-1 zeros between every element.
+## @deftypefn {Function File} {@var{y} =} upsample (@var{x}, @var{n})
+## @deftypefnx {Function File} {@var{y} =} upsample (@var{x}, @var{n}, @var{offset})
+## Upsample the signal, inserting n-1 zeros between every element.
+##
## If @var{x} is a matrix, upsample every column.
##
-## @deftypefnx {Function File} @var{y} = upsample(@var{x},@var{n},@var{offset})
-## Control the position of the inserted sample in the block of n zeros.
+## If @var{offset} is specified, control the position of the inserted sample in
+## the block of n zeros.
+## @seealso{decimate, downsample, interp, resample, upfirdn}
## @end deftypefn
-## @seealso{decimate, downsample, interp, resample, upfirdn}
-function y = upsample(x,n,phase)
- if nargin<2 || nargin>3, usage('upsample(x,n,[phase]'); end
- if nargin==2, phase = 0; end
+function y = upsample (x, n, phase = 0)
+ if nargin<2 || nargin>3, print_usage; end
if phase > n - 1
- warning("This is incompatible with Matlab (phase = 0:n-1). See\
+ warning("This is incompatible with Matlab (phase = 0:n-1). See ...
octave-forge signal package release notes for details." )
end
@@ -36,4 +37,3 @@
%!assert(upsample([1,2;5,6;9,10],2),[1,2;0,0;5,6;0,0;9,10;0,0]);
%!assert(upsample([2,4],2,1),[0,2,0,4]);
%!assert(upsample([3,4;7,8],2,1),[0,0;3,4;0,0;7,8]);
-
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-01-06 14:29:02
|
Revision: 9500
http://octave.svn.sourceforge.net/octave/?rev=9500&view=rev
Author: paramaniac
Date: 2012-01-06 14:28:55 +0000 (Fri, 06 Jan 2012)
Log Message:
-----------
signal: add x2y conversions
Added Paths:
-----------
trunk/octave-forge/main/signal/inst/ss2tf.m
trunk/octave-forge/main/signal/inst/ss2zp.m
trunk/octave-forge/main/signal/inst/tf2ss.m
trunk/octave-forge/main/signal/inst/tf2zp.m
trunk/octave-forge/main/signal/inst/zp2ss.m
trunk/octave-forge/main/signal/inst/zp2tf.m
Added: trunk/octave-forge/main/signal/inst/ss2tf.m
===================================================================
--- trunk/octave-forge/main/signal/inst/ss2tf.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/ss2tf.m 2012-01-06 14:28:55 UTC (rev 9500)
@@ -0,0 +1,66 @@
+## Copyright (C) 1996, 2000, 2004, 2005, 2007
+## Auburn University. All rights reserved.
+## Copyright (C) 2012 Lukas F. Reichlin
+##
+## This program is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{num}, @var{den}] =} ss2tf (@var{a}, @var{b}, @var{c}, @var{d})
+## Conversion from transfer function to state-space.
+## The state space system:
+## @iftex
+## @tex
+## $$ \dot x = Ax + Bu $$
+## $$ y = Cx + Du $$
+## @end tex
+## @end iftex
+## @ifinfo
+## @example
+## .
+## x = Ax + Bu
+## y = Cx + Du
+## @end example
+## @end ifinfo
+##
+## is converted to a transfer function:
+## @iftex
+## @tex
+## $$ G(s) = { { \rm num }(s) \over { \rm den }(s) } $$
+## @end tex
+## @end iftex
+## @ifinfo
+## @example
+##
+## num(s)
+## G(s)=-------
+## den(s)
+## @end example
+## @end ifinfo
+##
+## @end deftypefn
+
+## Author: R. Bruce Tenison <bte...@en...>
+## Created: June 24, 1994
+## a s hodel: modified to allow for pure gain blocks Aug 1996
+
+function [num, den] = ss2tf (varargin)
+
+ if (nargin == 0)
+ print_usage ();
+ endif
+
+ [num, den] = tfdata (ss (varargin{:}), "vector");
+
+endfunction
Added: trunk/octave-forge/main/signal/inst/ss2zp.m
===================================================================
--- trunk/octave-forge/main/signal/inst/ss2zp.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/ss2zp.m 2012-01-06 14:28:55 UTC (rev 9500)
@@ -0,0 +1,39 @@
+## Copyright (C) 1996, 2000, 2004, 2005, 2006, 2007
+## Auburn University. All rights reserved.
+## Copyright (C) 2012 Lukas F. Reichlin
+##
+## This program is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{pol}, @var{zer}, @var{k}] =} ss2zp (@var{a}, @var{b}, @var{c}, @var{d})
+## Converts a state space representation to a set of poles and zeros;
+## @var{k} is a gain associated with the zeros.
+##
+## @end deftypefn
+
+## Author: David Clem
+## Created: August 15, 1994
+## Hodel: changed order of output arguments to zer, pol, k. July 1996
+## a s hodel: added argument checking, allow for pure gain blocks aug 1996
+
+function [z, p, k] = ss2zp (varargin)
+
+ if (nargin == 0)
+ print_usage ();
+ endif
+
+ [z, p, k] = zpkdata (ss (varargin{:}), "vector");
+
+endfunction
Added: trunk/octave-forge/main/signal/inst/tf2ss.m
===================================================================
--- trunk/octave-forge/main/signal/inst/tf2ss.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/tf2ss.m 2012-01-06 14:28:55 UTC (rev 9500)
@@ -0,0 +1,67 @@
+## Copyright (C) 1996, 1998, 2000, 2002, 2004, 2005, 2007
+## Auburn University. All rights reserved.
+## Copyright (C) 2012 Lukas F. Reichlin
+##
+## This program is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{a}, @var{b}, @var{c}, @var{d}] =} tf2ss (@var{num}, @var{den})
+## Conversion from transfer function to state-space.
+## The state space system:
+## @iftex
+## @tex
+## $$ \dot x = Ax + Bu $$
+## $$ y = Cx + Du $$
+## @end tex
+## @end iftex
+## @ifinfo
+## @example
+## .
+## x = Ax + Bu
+## y = Cx + Du
+## @end example
+## @end ifinfo
+## is obtained from a transfer function:
+## @iftex
+## @tex
+## $$ G(s) = { { \rm num }(s) \over { \rm den }(s) } $$
+## @end tex
+## @end iftex
+## @ifinfo
+## @example
+## num(s)
+## G(s)=-------
+## den(s)
+## @end example
+## @end ifinfo
+##
+## The state space system matrices obtained from this function
+## will be in observable companion form as Wolovich's Observable
+## Structure Theorem is used.
+## @end deftypefn
+
+## Author: R. Bruce Tenison <bte...@en...>
+## Created: June 22, 1994
+## mod A S Hodel July, Aug 1995
+
+function [a, b, c, d, e] = tf2ss (varargin)
+
+ if (nargin == 0)
+ print_usage ();
+ endif
+
+ [a, b, c, d, e] = dssdata (tf (varargin{:}), []);
+
+endfunction
Added: trunk/octave-forge/main/signal/inst/tf2zp.m
===================================================================
--- trunk/octave-forge/main/signal/inst/tf2zp.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/tf2zp.m 2012-01-06 14:28:55 UTC (rev 9500)
@@ -0,0 +1,38 @@
+## Copyright (C) 1996, 1998, 2000, 2003, 2004, 2005, 2006, 2007
+## Auburn University. All rights reserved.
+## Copyright (C) 2012 Lukas F. Reichlin
+##
+## This program is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{zer}, @var{pol}, @var{k}] =} tf2zp (@var{num}, @var{den})
+## Converts transfer functions to poles-and-zero representations.
+##
+## Returns the zeros and poles of the system defined
+## by @var{num}/@var{den}.
+## @var{k} is a gain associated with the system zeros.
+## @end deftypefn
+
+## Author: A. S. Hodel <a.s...@en...>
+
+function [z, p, k] = tf2zp (varargin)
+
+ if (nargin == 0)
+ print_usage ();
+ endif
+
+ [z, p, k] = zpkdata (tf (varargin{:}), "vector");
+
+endfunction
Added: trunk/octave-forge/main/signal/inst/zp2ss.m
===================================================================
--- trunk/octave-forge/main/signal/inst/zp2ss.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/zp2ss.m 2012-01-06 14:28:55 UTC (rev 9500)
@@ -0,0 +1,69 @@
+## Copyright (C) 1996, 2000, 2002, 2003, 2004, 2005, 2007
+## Auburn University. All rights reserved.
+## Copyright (C) 2012 Lukas F. Reichlin
+##
+## This program is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{a}, @var{b}, @var{c}, @var{d}] =} zp2ss (@var{zer}, @var{pol}, @var{k})
+## Conversion from zero / pole to state space.
+##
+## @strong{Inputs}
+## @table @var
+## @item zer
+## @itemx pol
+## Vectors of (possibly) complex poles and zeros of a transfer
+## function. Complex values must come in conjugate pairs
+## (i.e., @math{x+jy} in @var{zer} means that @math{x-jy} is also in @var{zer}).
+## @item k
+## Real scalar (leading coefficient).
+## @end table
+##
+## @strong{Outputs}
+## @table @var
+## @item @var{a}
+## @itemx @var{b}
+## @itemx @var{c}
+## @itemx @var{d}
+## The state space system, in the form:
+## @iftex
+## @tex
+## $$ \dot x = Ax + Bu $$
+## $$ y = Cx + Du $$
+## @end tex
+## @end iftex
+## @ifinfo
+## @example
+## .
+## x = Ax + Bu
+## y = Cx + Du
+## @end example
+## @end ifinfo
+## @end table
+## @end deftypefn
+
+## Author: David Clem
+## Created: August 15, 1994
+
+function [a, b, c, d, e] = zp2ss (varargin)
+
+ if (nargin == 0)
+ print_usage ();
+ endif
+
+ [a, b, c, d, e] = dssdata (zpk (varargin{:}), []);
+
+endfunction
+
Added: trunk/octave-forge/main/signal/inst/zp2tf.m
===================================================================
--- trunk/octave-forge/main/signal/inst/zp2tf.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/zp2tf.m 2012-01-06 14:28:55 UTC (rev 9500)
@@ -0,0 +1,45 @@
+## Copyright (C) 1996, 1998, 2000, 2002, 2004, 2005, 2007
+## Auburn University. All rights reserved.
+## Copyright (C) 2012 Lukas F. Reichlin
+##
+## This program is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{num}, @var{den}] =} zp2tf (@var{zer}, @var{pol}, @var{k})
+## Converts zeros / poles to a transfer function.
+##
+## @strong{Inputs}
+## @table @var
+## @item zer
+## @itemx pol
+## Vectors of (possibly complex) poles and zeros of a transfer
+## function. Complex values must appear in conjugate pairs.
+## @item k
+## Real scalar (leading coefficient).
+## @end table
+## @end deftypefn
+
+## Author: A. S. Hodel <a.s...@en...>
+## (With help from students Ingram, McGowan.)
+
+function [num, den] = zp2tf (varargin)
+
+ if (nargin == 0)
+ print_usage ();
+ endif
+
+ [num, den] = tfdata (zpk (varargin{:}), "vector");
+
+endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <car...@us...> - 2012-01-18 14:03:07
|
Revision: 9531
http://octave.svn.sourceforge.net/octave/?rev=9531&view=rev
Author: carandraug
Date: 2012-01-18 14:02:54 +0000 (Wed, 18 Jan 2012)
Log Message:
-----------
signal:
* update GPL to v3 (or later)
* fix copyright notices
* use print_usage
* fix indentation
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/__ellip_ws.m
trunk/octave-forge/main/signal/inst/__ellip_ws_min.m
trunk/octave-forge/main/signal/inst/arburg.m
trunk/octave-forge/main/signal/inst/besself.m
trunk/octave-forge/main/signal/inst/bilinear.m
trunk/octave-forge/main/signal/inst/bitrevorder.m
trunk/octave-forge/main/signal/inst/blackmanharris.m
trunk/octave-forge/main/signal/inst/blackmannuttall.m
trunk/octave-forge/main/signal/inst/bohmanwin.m
trunk/octave-forge/main/signal/inst/boxcar.m
trunk/octave-forge/main/signal/inst/buffer.m
trunk/octave-forge/main/signal/inst/butter.m
trunk/octave-forge/main/signal/inst/buttord.m
trunk/octave-forge/main/signal/inst/cceps.m
trunk/octave-forge/main/signal/inst/cheb.m
trunk/octave-forge/main/signal/inst/cheb1ord.m
trunk/octave-forge/main/signal/inst/cheb2ord.m
trunk/octave-forge/main/signal/inst/chebwin.m
trunk/octave-forge/main/signal/inst/cheby1.m
trunk/octave-forge/main/signal/inst/cheby2.m
trunk/octave-forge/main/signal/inst/chirp.m
trunk/octave-forge/main/signal/inst/cmorwavf.m
trunk/octave-forge/main/signal/inst/cohere.m
trunk/octave-forge/main/signal/inst/convmtx.m
trunk/octave-forge/main/signal/inst/cplxreal.m
trunk/octave-forge/main/signal/inst/cpsd.m
trunk/octave-forge/main/signal/inst/csd.m
trunk/octave-forge/main/signal/inst/czt.m
trunk/octave-forge/main/signal/inst/dct.m
trunk/octave-forge/main/signal/inst/dct2.m
trunk/octave-forge/main/signal/inst/dctmtx.m
trunk/octave-forge/main/signal/inst/decimate.m
trunk/octave-forge/main/signal/inst/dftmtx.m
trunk/octave-forge/main/signal/inst/diric.m
trunk/octave-forge/main/signal/inst/downsample.m
trunk/octave-forge/main/signal/inst/dst.m
trunk/octave-forge/main/signal/inst/dwt.m
trunk/octave-forge/main/signal/inst/ellip.m
trunk/octave-forge/main/signal/inst/ellipord.m
trunk/octave-forge/main/signal/inst/fht.m
trunk/octave-forge/main/signal/inst/filtfilt.m
trunk/octave-forge/main/signal/inst/filtic.m
trunk/octave-forge/main/signal/inst/fir1.m
trunk/octave-forge/main/signal/inst/fir2.m
Modified: trunk/octave-forge/main/signal/inst/__ellip_ws.m
===================================================================
--- trunk/octave-forge/main/signal/inst/__ellip_ws.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/__ellip_ws.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
Modified: trunk/octave-forge/main/signal/inst/__ellip_ws_min.m
===================================================================
--- trunk/octave-forge/main/signal/inst/__ellip_ws_min.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/__ellip_ws_min.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -22,14 +22,12 @@
##
## - Serra, Celso Penteado, Teoria e Projeto de Filtros, Campinas: CARTGRAF,
## 1983.
-## Author: Paulo Neis <p_...@ya...>
function err=__ellip_ws_min(kl, x)
-#
-int=ellipke([kl; 1-kl]);
-ql=int(1);
-q=int(2);
-err=abs((ql/q)-x);
+ int=ellipke([kl; 1-kl]);
+ ql=int(1);
+ q=int(2);
+ err=abs((ql/q)-x);
endfunction
Modified: trunk/octave-forge/main/signal/inst/arburg.m
===================================================================
--- trunk/octave-forge/main/signal/inst/arburg.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/arburg.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License
-%% as published by the Free Software Foundation; either version 2,
+%% as published by the Free Software Foundation; either version 3,
%% or (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
Modified: trunk/octave-forge/main/signal/inst/besself.m
===================================================================
--- trunk/octave-forge/main/signal/inst/besself.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/besself.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -1,5 +1,6 @@
-## Copyright (C) 2009 Thomas Sailer
## Copyright (C) 1999 Paul Kienzle <pki...@us...>
+## Copyright (C) 2003 Doug Stewart <da...@sy...>
+## Copyright (C) 2009 Thomas Sailer <t.s...@al...>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
@@ -44,9 +45,6 @@
## Proakis & Manolakis (1992). Digital Signal Processing. New York:
## Macmillan Publishing Company.
-## Author: Paul Kienzle <pki...@us...>
-## Modified by: Doug Stewart <da...@sy...> Feb, 2003
-
function [a, b, c, d] = besself (n, W, varargin)
if (nargin>4 || nargin<2) || (nargout>4 || nargout<2)
Modified: trunk/octave-forge/main/signal/inst/bilinear.m
===================================================================
--- trunk/octave-forge/main/signal/inst/bilinear.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/bilinear.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -70,8 +70,6 @@
## Proakis & Manolakis (1992). Digital Signal Processing. New York:
## Macmillan Publishing Company.
-## Author: Paul Kienzle <pki...@us...>
-
function [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T)
if nargin==3
Modified: trunk/octave-forge/main/signal/inst/bitrevorder.m
===================================================================
--- trunk/octave-forge/main/signal/inst/bitrevorder.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/bitrevorder.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -1,4 +1,4 @@
-## Copyright (C) 2007 Sylvain Pelissier <syl...@gm...>
+## Copyright (C) 2007 Sylvain Pelissier <syl...@gm...>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
@@ -15,7 +15,7 @@
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{y} @var{i}] =} bitrevorder(@var{x})
-## Reorder x in the bit reversed order
+## Reorder x in the bit reversed order
## @seealso{fft,ifft}
## @end deftypefn
@@ -33,4 +33,4 @@
new_ind = bi2de(fliplr(de2bi(old_ind)));
i = new_ind + 1;
- y(old_ind+1) = x(i);
\ No newline at end of file
+ y(old_ind+1) = x(i);
Modified: trunk/octave-forge/main/signal/inst/blackmanharris.m
===================================================================
--- trunk/octave-forge/main/signal/inst/blackmanharris.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/blackmanharris.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -1,4 +1,4 @@
-## Copyright (C) 2007 Sylvain Pelissier <syl...@gm...>
+## Copyright (C) 2007 Sylvain Pelissier <syl...@gm...>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
@@ -15,22 +15,22 @@
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{w}] =} blackmanharris(@var{L})
-## Compute the Blackman-Harris window.
+## Compute the Blackman-Harris window.
## @seealso{rectwin, bartlett}
## @end deftypefn
-function [w] = blackmanharris(L)
- if (nargin < 1); usage('blackmanharris(x)'); end
- if(! isscalar(L))
- error("L must be a number");
- endif
-
- N = L-1;
- a0 = 0.35875;
- a1 = 0.48829;
- a2 = 0.14128;
- a3 = 0.01168;
- n = -ceil(N/2):N/2;
- w = a0 + a1.*cos(2.*pi.*n./N) + a2.*cos(4.*pi.*n./N) + a3.*cos(6.*pi.*n./N);
-endfunction;
-
\ No newline at end of file
+function [w] = blackmanharris (L)
+ if (nargin < 1)
+ print_usage;
+ elseif(! isscalar(L))
+ error("L must be a number");
+ endif
+
+ N = L-1;
+ a0 = 0.35875;
+ a1 = 0.48829;
+ a2 = 0.14128;
+ a3 = 0.01168;
+ n = -ceil(N/2):N/2;
+ w = a0 + a1.*cos(2.*pi.*n./N) + a2.*cos(4.*pi.*n./N) + a3.*cos(6.*pi.*n./N);
+endfunction
Modified: trunk/octave-forge/main/signal/inst/blackmannuttall.m
===================================================================
--- trunk/octave-forge/main/signal/inst/blackmannuttall.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/blackmannuttall.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -1,4 +1,4 @@
-## Copyright (C) 2007 Muthiah Annamalai <mut...@ut...>
+## Copyright (C) 2007 Muthiah Annamalai <mut...@ut...>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
@@ -15,22 +15,23 @@
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{w}] =} blackmannuttall(@var{L})
-## Compute the Blackman-Nuttall window.
+## Compute the Blackman-Nuttall window.
## @seealso{nuttallwin, kaiser}
## @end deftypefn
function [w] = blackmannuttall(L)
- if (nargin < 1); usage('blackmannuttall(L)'); end
- if(! isscalar(L))
- error("L must be a number");
- endif
-
- N = L-1;
- a0 = 0.3635819;
- a1 = 0.4891775;
- a2 = 0.1365995;
- a3 = 0.0106411;
- n = 0:N;
- w = a0 - a1.*cos(2.*pi.*n./N) + a2.*cos(4.*pi.*n./N) - a3.*cos(6.*pi.*n./N);
- w = w.';
-endfunction;
+ if (nargin < 1)
+ print_usage;
+ elseif (! isscalar(L))
+ error("L must be a number");
+ endif
+
+ N = L-1;
+ a0 = 0.3635819;
+ a1 = 0.4891775;
+ a2 = 0.1365995;
+ a3 = 0.0106411;
+ n = 0:N;
+ w = a0 - a1.*cos(2.*pi.*n./N) + a2.*cos(4.*pi.*n./N) - a3.*cos(6.*pi.*n./N);
+ w = w.';
+endfunction
Modified: trunk/octave-forge/main/signal/inst/bohmanwin.m
===================================================================
--- trunk/octave-forge/main/signal/inst/bohmanwin.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/bohmanwin.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -1,4 +1,4 @@
-## Copyright (C) 2007 Sylvain Pelissier <syl...@gm...>
+## Copyright (C) 2007 Sylvain Pelissier <syl...@gm...>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
@@ -15,38 +15,37 @@
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{w}] =} bohmanwin(@var{L})
-## Compute the Bohman window of lenght L.
+## Compute the Bohman window of lenght L.
## @seealso{rectwin, bartlett}
## @end deftypefn
function [w] = bohmanwin(L)
- if (nargin < 1); usage('bohmanwin(x)'); end
- if(! isscalar(L))
- error("L must be a number");
- endif
-
- if(L < 0)
- error('L must be positive');
- end
-
- if(L ~= floor(L))
- L = round(L);
- warning('L rounded to the nearest integer.');
- end
-
- if(L == 0)
- w = [];
-
- elseif(L == 1)
- w = 1;
-
- else
- N = L-1;
- n = -N/2:N/2;
-
- w = (1-2.*abs(n)./N).*cos(2.*pi.*abs(n)./N) + (1./pi).*sin(2.*pi.*abs(n)./N);
- w(1) = 0;
- w(length(w))=0;
- w = w';
- end
-endfunction;
\ No newline at end of file
+ if (nargin < 1)
+ print_usage
+ elseif(! isscalar(L))
+ error("L must be a number");
+ elseif(L < 0)
+ error('L must be positive');
+ end
+
+ if(L ~= floor(L))
+ L = round(L);
+ warning('L rounded to the nearest integer.');
+ end
+
+ if(L == 0)
+ w = [];
+
+ elseif(L == 1)
+ w = 1;
+
+ else
+ N = L-1;
+ n = -N/2:N/2;
+
+ w = (1-2.*abs(n)./N).*cos(2.*pi.*abs(n)./N) + (1./pi).*sin(2.*pi.*abs(n)./N);
+ w(1) = 0;
+ w(length(w))=0;
+ w = w';
+ end
+endfunction
Modified: trunk/octave-forge/main/signal/inst/boxcar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/boxcar.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/boxcar.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -18,15 +18,13 @@
## Returns the filter coefficients of a rectangular window of length n.
function w = boxcar (n)
-
+
if (nargin != 1)
- usage ("w = boxcar(n)");
- endif
-
- if !isscalar(n) || n != floor(n) || n <= 0
+ print_usage;
+ elseif !isscalar(n) || n != floor(n) || n <= 0
error ("boxcar: n must be an integer > 0");
endif
w = ones(n, 1);
-
+
endfunction
Modified: trunk/octave-forge/main/signal/inst/buffer.m
===================================================================
--- trunk/octave-forge/main/signal/inst/buffer.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/buffer.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
Modified: trunk/octave-forge/main/signal/inst/butter.m
===================================================================
--- trunk/octave-forge/main/signal/inst/butter.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/butter.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -1,8 +1,10 @@
## Copyright (C) 1999 Paul Kienzle <pki...@us...>
+## Copyright (C) 2003 Doug Stewart <da...@sy...>
+## Copyright (C) 2011 Alexander Klein <ale...@ma...>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -43,14 +45,10 @@
## Proakis & Manolakis (1992). Digital Signal Processing. New York:
## Macmillan Publishing Company.
-## Author: Paul Kienzle <pki...@us...>
-## Modified by: Doug Stewart <da...@sy...> Feb, 2003
-## Tests and demos added: Alexander Klein <ale...@ma...>, Sep 2011
-
function [a, b, c, d] = butter (n, W, varargin)
if (nargin>4 || nargin<2) || (nargout>4 || nargout<2)
- usage ("[b, a] or [z, p, g] or [a,b,c,d] = butter (n, W [, 'ftype'][,'s'])");
+ print_usage;
end
## interpret the input parameters
@@ -59,7 +57,7 @@
end
stop = 0;
- digital = 1;
+ digital = 1;
for i=1:length(varargin)
switch varargin{i}
case 's', digital = 0;
Modified: trunk/octave-forge/main/signal/inst/buttord.m
===================================================================
--- trunk/octave-forge/main/signal/inst/buttord.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/buttord.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -39,15 +39,12 @@
function [n, Wc] = buttord(Wp, Ws, Rp, Rs)
if nargin != 4
- usage("[n, Wn] = buttord(Wp, Ws, Rp, Rs)");
- end
- if length(Wp) != length(Ws)
+ print_usage;
+ elseif length(Wp) != length(Ws)
error("buttord: Wp and Ws must have the same length");
- end
- if length(Wp) != 1 && length(Wp) != 2
+ elseif length(Wp) != 1 && length(Wp) != 2
error("buttord: Wp,Ws must have length 1 or 2");
- end
- if length(Wp) == 2 && (all(Wp>Ws) || all(Ws>Wp) || diff(Wp)<=0 || diff(Ws)<=0)
+ elseif length(Wp) == 2 && (all(Wp>Ws) || all(Ws>Wp) || diff(Wp)<=0 || diff(Ws)<=0)
error("buttord: Wp(1)<Ws(1)<Ws(2)<Wp(2) or Ws(1)<Wp(1)<Wp(2)<Ws(2)");
end
@@ -83,4 +80,3 @@
Wc(stop) = 1-Wc(stop);
endfunction
-
Modified: trunk/octave-forge/main/signal/inst/cceps.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cceps.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/cceps.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -22,13 +22,13 @@
## Author: Andreas Weingessel <And...@ci...>
## Apr 1, 1994
## Last modifified by AW on Nov 8, 1994
-
+
function cep = cceps (x, c)
if (nargin == 1)
c = 0;
elseif (nargin != 2)
- error ("usage: cceps (x [, correct])");
+ print_usage;
endif
[nr, nc] = size (x);
@@ -57,7 +57,7 @@
if (cor)
F = fft (x(1:nr-1))
if (min (abs (F)) == 0)
- error (bad_signal_message);
+ error (bad_signal_message);
endif
endif
endif
Modified: trunk/octave-forge/main/signal/inst/cheb.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cheb.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/cheb.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or (at
+## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## This program is distributed in the hope that it will be useful, but
@@ -25,22 +25,17 @@
## If x is a vector, the output is a vector of the same size, where each
## element is calculated as y(i) = Tn(x(i)).
-## Author: André Carezia <aca...@uo...>
-## Description: Value of the Chebyshev polynomials
-
function T = cheb (n, x)
if (nargin != 2)
- usage ("cheb (n, x)");
- endif
-
- if !(isscalar (n) && (n == round(n)) && (n >= 0))
+ print_usage;
+ elseif !(isscalar (n) && (n == round(n)) && (n >= 0))
error ("cheb: n has to be a positive integer");
endif
if (max(size(x)) == 0)
T = [];
endif
- # avoid resizing latencies
+ # avoid resizing latencies
T = zeros(size(x));
ind = abs (x) <= 1;
if (max(size(ind)))
@@ -53,5 +48,4 @@
endif
T = real(T);
-end
-
+endfunction
Modified: trunk/octave-forge/main/signal/inst/cheb1ord.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cheb1ord.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/cheb1ord.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -1,8 +1,9 @@
## Copyright (C) 2000 Paul Kienzle <pki...@us...>
+## Copyright (C) 2000 Laurent S. Mazet
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -12,8 +13,6 @@
##
## You should have received a copy of the GNU General Public License
## along with this program; If not, see <http://www.gnu.org/licenses/>.
-##
-## Completed by: Laurent S. Mazet
## Compute chebyshev type I filter order and cutoff for the desired response
## characteristics. Rp is the allowable decibels of ripple in the pass
@@ -34,13 +33,13 @@
function [n, Wc] = cheb1ord(Wp, Ws, Rp, Rs)
if nargin != 4
- usage("[n, Wn] = cheb1ord(Wp, Ws, Rp, Rs)");
+ print_usage;
elseif length(Wp) != length(Ws)
error("cheb1ord: Wp and Ws must have the same length");
elseif length(Wp) != 1 && length(Wp) != 2
error("cheb1ord: Wp,Ws must have length 1 or 2");
elseif length(Wp) == 2 && ...
- (all(Wp>Ws) || all(Ws>Wp) || diff(Wp)<=0 || diff(Ws)<=0)
+ (all(Wp>Ws) || all(Ws>Wp) || diff(Wp)<=0 || diff(Ws)<=0)
error("cheb1ord: Wp(1)<Ws(1)<Ws(2)<Wp(2) or Ws(1)<Wp(1)<Wp(2)<Ws(2)");
end
Modified: trunk/octave-forge/main/signal/inst/cheb2ord.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cheb2ord.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/cheb2ord.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
Modified: trunk/octave-forge/main/signal/inst/chebwin.m
===================================================================
--- trunk/octave-forge/main/signal/inst/chebwin.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/chebwin.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or (at
+## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## This program is distributed in the hope that it will be useful, but
@@ -48,21 +48,13 @@
##
## See also: kaiser
-## $Id$
-##
-## Author: André Carezia <aca...@uo...>
-## Description: Coefficients of the Dolph-Chebyshev window
-
function w = chebwin (n, at)
if (nargin != 2)
- usage ("chebwin (n, at)");
- endif
-
- if !(isscalar (n) && (n == round(n)) && (n > 0))
+ print_usage;
+ elseif !(isscalar (n) && (n == round(n)) && (n > 0))
error ("chebwin: n has to be a positive integer");
- endif
- if !(isscalar (at) && (at == real (at)))
+ elseif !(isscalar (at) && (at == real (at)))
error ("chebwin: at has to be a real scalar");
endif
@@ -94,5 +86,4 @@
endif
w = w ./ max (w (:));
-end
-
+endfunction
Modified: trunk/octave-forge/main/signal/inst/cheby1.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cheby1.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/cheby1.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -1,8 +1,9 @@
## Copyright (C) 1999 Paul Kienzle <pki...@us...>
+## Copyright (C) 2003 Doug Stewart <da...@sy...>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -42,13 +43,10 @@
## Parks & Burrus (1987). Digital Filter Design. New York:
## John Wiley & Sons, Inc.
-## Author: Paul Kienzle <pki...@us...>
-## Modified: Doug Stewart Feb. 2003
-
function [a,b,c,d] = cheby1(n, Rp, W, varargin)
if (nargin>5 || nargin<3) || (nargout>4 || nargout<2)
- usage ("[b, a] or [z, p, g] or [a,b,c,d]= cheby1 (n, Rp, W, [, 'ftype'][,'s'])");
+ print_usage;
endif
## interpret the input parameters
Modified: trunk/octave-forge/main/signal/inst/cheby2.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cheby2.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/cheby2.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -1,8 +1,9 @@
## Copyright (C) 1999 Paul Kienzle <pki...@us...>
+## Copyright (C) 2003 Doug Stewart <da...@sy...>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -42,13 +43,10 @@
## Parks & Burrus (1987). Digital Filter Design. New York:
## John Wiley & Sons, Inc.
-## Author: Paul Kienzle <pki...@us...>
-## Modified: Doug Stewart Feb. 2003
-
function [a,b,c,d] = cheby2(n, Rs, W, varargin)
if (nargin>5 || nargin<3) || (nargout>4 || nargout<2)
- usage ("[b, a] or [z, p, g] or [a,b,c,d]= cheby2 (n, Rs, W, [, 'ftype'][,'s'])");
+ print_usage;
end
## interpret the input parameters
Modified: trunk/octave-forge/main/signal/inst/chirp.m
===================================================================
--- trunk/octave-forge/main/signal/inst/chirp.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/chirp.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -36,13 +36,10 @@
## If you want a different sweep shape f(t), use the following:
## y = cos(2*pi*integral(f(t)) + 2*pi*f0*t + phase);
-## 2001-08-31 Paul Kienzle <pki...@us...>
-## * Fix documentation for quadratic case
-
function y = chirp(t, f0, t1, f1, form, phase)
if nargin < 1 || nargin > 6
- usage("y = chirp(t [, f0 [, t1 [, f1 [, form [, phase]]]]])");
+ print_usage;
endif
if nargin < 2, f0 = []; endif
if nargin < 3, t1 = []; endif
Modified: trunk/octave-forge/main/signal/inst/cmorwavf.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cmorwavf.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/cmorwavf.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -1,8 +1,8 @@
-## Copyright (C) 2007 Sylvain Pelissier <syl...@gm...>
+## Copyright (C) 2007 Sylvain Pelissier <syl...@gm...>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -16,16 +16,15 @@
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{psi,x}] =} cmorwavf (@var{lb,ub,n,fb,fc})
-## Compute the Complex Morlet wavelet.
+## Compute the Complex Morlet wavelet.
## @end deftypefn
function [psi,x] = cmorwavf (lb,ub,n,fb,fc)
- if (nargin ~= 5); usage('[psi,x] = cmorwavf(lb,ub,n,fb,fc)'); end
-
- if (n <= 0 || floor(n) ~= n)
- error("n must be an integer strictly positive");
- endif
- x = linspace(lb,ub,n);
- psi =((pi*fb)^(-0.5))*exp(2*i*pi*fc.*x).*exp(-x.^2/fb);
+ if (nargin ~= 5)
+ print_usage;
+ elseif (n <= 0 || floor(n) ~= n)
+ error("n must be an integer strictly positive");
+ endif
+ x = linspace(lb,ub,n);
+ psi =((pi*fb)^(-0.5))*exp(2*i*pi*fc.*x).*exp(-x.^2/fb);
endfunction
-
Modified: trunk/octave-forge/main/signal/inst/cohere.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cohere.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/cohere.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License
-%% as published by the Free Software Foundation; either version 2,
+%% as published by the Free Software Foundation; either version 3,
%% or (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
@@ -22,9 +22,7 @@
%% See "help pwelch" for description of arguments, hints and references
%% --- especially hint (7) for Matlab R11 defaults.
%%
-%%
-
function [varargout] = cohere(varargin)
%%
if ( nargin<2 )
Modified: trunk/octave-forge/main/signal/inst/convmtx.m
===================================================================
--- trunk/octave-forge/main/signal/inst/convmtx.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/convmtx.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -37,7 +37,7 @@
function b = convmtx (a, n)
if (nargin != 2)
- usage ("convmtx (a, n)");
+ print_usage;
endif
[r, c] = size(a);
Modified: trunk/octave-forge/main/signal/inst/cplxreal.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cplxreal.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/cplxreal.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
%%
%% This program is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
-%% the Free Software Foundation; either version 2 of the License, or
+%% the Free Software Foundation; either version 3 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
Modified: trunk/octave-forge/main/signal/inst/cpsd.m
===================================================================
--- trunk/octave-forge/main/signal/inst/cpsd.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/cpsd.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License
-%% as published by the Free Software Foundation; either version 2,
+%% as published by the Free Software Foundation; either version 3,
%% or (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
@@ -19,9 +19,7 @@
%% Estimate cross power spectrum of data "x" and "y" by the Welch (1967)
%% periodogram/FFT method.
%% See "help pwelch" for description of arguments, hints and references
-%%
-
function [varargout] = cpsd(varargin)
%%
%% Check fixed argument
Modified: trunk/octave-forge/main/signal/inst/csd.m
===================================================================
--- trunk/octave-forge/main/signal/inst/csd.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/csd.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License
-%% as published by the Free Software Foundation; either version 2,
+%% as published by the Free Software Foundation; either version 3,
%% or (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
@@ -21,9 +21,7 @@
%% See "help pwelch" for description of arguments, hints and references
%% --- especially hint (7) for Matlab R11 defaults.
%%
-%%
-
function [varargout] = csd(varargin)
%%
%% Check fixed argument
Modified: trunk/octave-forge/main/signal/inst/czt.m
===================================================================
--- trunk/octave-forge/main/signal/inst/czt.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/czt.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -44,7 +44,7 @@
## multiply gg by M-elements of chirp and call it done
function y = czt(x, m, w, a)
- if nargin < 1 || nargin > 4, usage("y=czt(x, m, w, a)"); endif
+ if nargin < 1 || nargin > 4, print_usage; endif
[row, col] = size(x);
if row == 1, x = x(:); col = 1; endif
Modified: trunk/octave-forge/main/signal/inst/dct.m
===================================================================
--- trunk/octave-forge/main/signal/inst/dct.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/dct.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -54,13 +54,11 @@
##
## Scaling the result by w(k)/2 will give us the desired output.
-## Author: Paul Kienzle
-## 2001-02-08
-## * initial release
+
function y = dct (x, n)
if (nargin < 1 || nargin > 2)
- usage ("y = dct(x [, n])");
+ print_usage;
endif
realx = isreal(x);
Modified: trunk/octave-forge/main/signal/inst/dct2.m
===================================================================
--- trunk/octave-forge/main/signal/inst/dct2.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/dct2.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -20,14 +20,10 @@
## Computes the 2-D DCT of x after padding or trimming rows to m and
## columns to n.
-## Author: Paul Kienzle
-## 2001-02-08
-## * initial revision
-
function y = dct2 (x, m, n)
if (nargin < 1 || nargin > 3)
- usage("dct (x) or dct (x, m, n) or dct (x, [m n])");
+ print_usage;
endif
if nargin == 1
Modified: trunk/octave-forge/main/signal/inst/dctmtx.m
===================================================================
--- trunk/octave-forge/main/signal/inst/dctmtx.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/dctmtx.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -14,7 +14,7 @@
## along with this program; If not, see <http://www.gnu.org/licenses/>.
## T = dctmtx (n)
-## Return the DCT transformation matrix of size n x n.
+## Return the DCT transformation matrix of size n x n.
##
## If A is an n x n matrix, then the following are true:
## T*A == dct(A), T'*A == idct(A)
@@ -29,13 +29,9 @@
##
## See also: dct, idct, dct2, idct2
-## Author: Paul Kienzle <pki...@us...>
-## 2001-02-08
-## * initial release
-
function T = dctmtx(n)
if nargin != 1
- usage("T = dctmtx(n)")
+ print_usage;
endif
if n > 1
Modified: trunk/octave-forge/main/signal/inst/decimate.m
===================================================================
--- trunk/octave-forge/main/signal/inst/decimate.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/decimate.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -35,10 +35,11 @@
function y = decimate(x, q, n, ftype)
- if nargin < 1 || nargin > 4,
- usage("y=decimate(x, q [, n] [, ftype])");
+ if nargin < 1 || nargin > 4
+ print_usage;
+ elseif q != fix(q)
+ error("decimate only works with integer q.");
endif
- if q != fix(q), error("decimate only works with integer q."); endif
if nargin<3
ftype='iir';
Modified: trunk/octave-forge/main/signal/inst/dftmtx.m
===================================================================
--- trunk/octave-forge/main/signal/inst/dftmtx.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/dftmtx.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -26,13 +26,11 @@
function d = dftmtx(n)
if (nargin != 1)
- error ("usage: d = dftmtx (n)");
- endif
-
- if (!isscalar(n))
+ print_usage;
+ elseif (!isscalar(n))
error ("dftmtx: argument must be scalar");
endif
-
+
d = fft(eye(n));
endfunction
Modified: trunk/octave-forge/main/signal/inst/diric.m
===================================================================
--- trunk/octave-forge/main/signal/inst/diric.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/diric.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -1,8 +1,8 @@
-## Copyright (C) 2007 Sylvain Pelissier <syl...@gm...>
+## Copyright (C) 2007 Sylvain Pelissier <syl...@gm...>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -15,18 +15,18 @@
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{y}] =} diric(@var{x},@var{n})
-## Compute the dirichlet function.
+## Compute the dirichlet function.
## @seealso{sinc, gauspuls, sawtooth}
## @end deftypefn
function [y] = diric(x,n)
- if (nargin < 2); usage('diric(x,n)'); end
-
- if (n <= 0 || floor(n) ~= n)
- error("n must be an integer strictly positive");
- endif
-
- y = sin(n.*x./2)./(n.*sin(x./2));
- y(mod(x,2*pi)==0) = (-1).^((n-1).*x(mod(x,2*pi)==0)./(2.*pi));
+ if (nargin < 2)
+ print_usage;
+ elseif (n <= 0 || floor(n) ~= n)
+ error("n must be an integer strictly positive");
+ endif
-endfunction;
\ No newline at end of file
+ y = sin(n.*x./2)./(n.*sin(x./2));
+ y(mod(x,2*pi)==0) = (-1).^((n-1).*x(mod(x,2*pi)==0)./(2.*pi));
+
+endfunction
Modified: trunk/octave-forge/main/signal/inst/downsample.m
===================================================================
--- trunk/octave-forge/main/signal/inst/downsample.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/downsample.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -1,4 +1,4 @@
-## Copyright (C) 2007 Paul Kienzle <pki...@us...>
+## Author: Paul Kienzle <pki...@us...> 2007
## This function is public domain
## -*- texinfo -*-
Modified: trunk/octave-forge/main/signal/inst/dst.m
===================================================================
--- trunk/octave-forge/main/signal/inst/dst.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/dst.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -1,4 +1,4 @@
-## Copyright (C) 2006 Paul Kienzle <pki...@us...>
+## Author: Paul Kienzle <pki...@us...> 2006
## This function is public domain
## -*- texinfo -*-
Modified: trunk/octave-forge/main/signal/inst/dwt.m
===================================================================
--- trunk/octave-forge/main/signal/inst/dwt.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/dwt.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -1,4 +1,4 @@
-## Copyright (C) 2008 Sylvain Pelissier <syl...@gm...>
+## Copyright (C) 2008 Sylvain Pelissier <syl...@gm...>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
@@ -15,19 +15,19 @@
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{ca} @var{cd}] =} dwt(@var{x,lo_d,hi_d})
-## Comupte de discrete wavelet transform of x with one level.
+## Comupte de discrete wavelet transform of x with one level.
## @end deftypefn
+function [ca cd] = dwt(x,lo_d,hi_d)
+ if (nargin < 3|| nargin > 3)
+ print_usage;
+ elseif(~isvector(x) || ~isvector(lo_d) || ~isvector(hi_d))
+ error('x, hi_d and lo_d must be vectors');
+ end
-function [ca cd] = dwt(x,lo_d,hi_d)
- if (nargin < 3|| nargin > 3); usage("[ca cd] = dwt(x,lo_d,hi_d)"); end
-
- if(~isvector(x) || ~isvector(lo_d) || ~isvector(hi_d))
- error('x, hi_d and lo_d must be vectors');
- end
-
- h = filter(hi_d,1,x);
- g = filter(lo_d,1,x);
-
- cd = downsample(h,2);
- ca = downsample(g,2);
\ No newline at end of file
+ h = filter(hi_d,1,x);
+ g = filter(lo_d,1,x);
+
+ cd = downsample(h,2);
+ ca = downsample(g,2);
+endfunction
Modified: trunk/octave-forge/main/signal/inst/ellip.m
===================================================================
--- trunk/octave-forge/main/signal/inst/ellip.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/ellip.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -1,8 +1,9 @@
## Copyright (C) 2001 Paulo Neis <p_...@ya...>
+## Copyright (C) 2003 Doug Stewart <da...@sy...>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -12,7 +13,6 @@
##
## You should have received a copy of the GNU General Public License
## along with this program; If not, see <http://www.gnu.org/licenses/>.
-##
## N-ellip 0.2.1
##usage: [Zz, Zp, Zg] = ellip(n, Rp, Rs, Wp, stype,'s')
@@ -47,15 +47,12 @@
## - Parente Ribeiro, E., Notas de aula da disciplina TE498 - Processamento
## Digital de Sinais, UFPR, 2001/2002.
## - Kienzle, Paul, functions from Octave-Forge, 1999 (http://octave.sf.net).
-##
-## Author: Paulo Neis <p_...@ya...>
-## Modified: Doug Stewart Feb. 2003
function [a,b,c,d] = ellip(n, Rp, Rs, W, varargin)
if (nargin>6 || nargin<4) || (nargout>4 || nargout<2)
- usage ("[b, a] or [z, p, g] or [a,b,c,d] = ellip (n, Rp, Rs, Wp, [, 'ftype'][,'s'])");
+ print_usage;
endif
## interpret the input parameters
Modified: trunk/octave-forge/main/signal/inst/ellipord.m
===================================================================
--- trunk/octave-forge/main/signal/inst/ellipord.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/ellipord.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -12,7 +12,6 @@
##
## You should have received a copy of the GNU General Public License
## along with this program; If not, see <http://www.gnu.org/licenses/>.
-##
## usage: [n,wp] = ellipord(wp,ws, rp,rs)
##
@@ -27,70 +26,65 @@
## - Lamar, Marcus Vinicius, Notas de aula da disciplina TE 456 - Circuitos
## Analogicos II, UFPR, 2001/2002.
-
-
function [n, Wp] = ellipord(Wp, Ws, Rp, Rs)
-##
-[rp, cp]=size(Wp);
-[rs, cs]=size(Ws);
-if ( !(length(Wp)<=2 && (rp==1 || cp==1) && length(Ws)<=2 && (rs==1 || cs==1)) )
- error ("ellipord: frequency must be given as w0 or [w0, w1]");
-elseif (!all(Wp >= 0 & Wp <= 1 & Ws >= 0 & Ws <= 1)) #####
- error ("ellipord: critical frequencies must be in (0 1)");
-elseif (!(length(Wp)==1 || length(Wp) == 2 & length(Ws)==1 || length(Ws) == 2))
- error ("ellipord: only one filter band allowed");
-elseif (length(Wp)==2 && !(Wp(1) < Wp(2)))
- error ("ellipord: first band edge must be smaller than second");
-elseif (length(Ws)==2 && !(length(Wp)==2))
- error ("ellipord: you must specify band pass borders.");
-elseif (length(Wp)==2 && length(Ws)==2 && !(Ws(1) < Wp(1) && Ws(2) > Wp(2)) )
- error ("ellipord: ( Wp(1), Wp(2) ) must be inside of interval ( Ws(1), Ws(2) )");
-elseif (length(Wp)==2 && length(Ws)==1 && !(Ws < Wp(1) || Ws > Wp(2)) )
- error ("ellipord: Ws must be out of interval ( Wp(1), Wp(2) )");
-endif
+ [rp, cp]=size(Wp);
+ [rs, cs]=size(Ws);
+ if ( !(length(Wp)<=2 && (rp==1 || cp==1) && length(Ws)<=2 && (rs==1 || cs==1)) )
+ error ("ellipord: frequency must be given as w0 or [w0, w1]");
+ elseif (!all(Wp >= 0 & Wp <= 1 & Ws >= 0 & Ws <= 1)) #####
+ error ("ellipord: critical frequencies must be in (0 1)");
+ elseif (!(length(Wp)==1 || length(Wp) == 2 & length(Ws)==1 || length(Ws) == 2))
+ error ("ellipord: only one filter band allowed");
+ elseif (length(Wp)==2 && !(Wp(1) < Wp(2)))
+ error ("ellipord: first band edge must be smaller than second");
+ elseif (length(Ws)==2 && !(length(Wp)==2))
+ error ("ellipord: you must specify band pass borders.");
+ elseif (length(Wp)==2 && length(Ws)==2 && !(Ws(1) < Wp(1) && Ws(2) > Wp(2)) )
+ error ("ellipord: ( Wp(1), Wp(2) ) must be inside of interval ( Ws(1), Ws(2) )");
+ elseif (length(Wp)==2 && length(Ws)==1 && !(Ws < Wp(1) || Ws > Wp(2)) )
+ error ("ellipord: Ws must be out of interval ( Wp(1), Wp(2) )");
+ endif
+ # sampling frequency of 2 Hz
+ T = 2;
+ Wpw = tan(pi.*Wp./T); # prewarp
+ Wsw = tan(pi.*Ws./T); # prewarp
+ ##pass/stop band to low pass filter transform:
+ if (length(Wpw)==2 && length(Wsw)==2)
+ wp=1;
+ w02 = Wpw(1) * Wpw(2); # Central frequency of stop/pass band (square)
+ w3 = w02/Wsw(2);
+ w4 = w02/Wsw(1);
+ if (w3 > Wsw(1))
+ ws = (Wsw(2)-w3)/(Wpw(2)-Wpw(1));
+ elseif (w4 < Wsw(2))
+ ws = (w4-Wsw(1))/(Wpw(2)-Wpw(1));
+ else
+ ws = (Wsw(2)-Wsw(1))/(Wpw(2)-Wpw(1));
+ endif
+ elseif (length(Wpw)==2 && length(Wsw)==1)
+ wp=1;
+ w02 = Wpw(1) * Wpw(2);
+ if (Wsw > Wpw(2))
+ w3 = w02/Wsw;
+ ws = (Wsw - w3)/(Wpw(2) - Wpw(1));
+ else
+ w4 = w02/Wsw;
+ ws = (w4 - Wsw)/(Wpw(2) - Wpw(1));
+ endif
+ else
+ wp = Wpw;
+ ws = Wsw;
+ endif
-# sampling frequency of 2 Hz
-T = 2;
+ k=wp/ws;
+ k1=sqrt(1-k^2);
+ q0=(1/2)*((1-sqrt(k1))/(1+sqrt(k1)));
+ q= q0 + 2*q0^5 + 15*q0^9 + 150*q0^13; %(....)
+ D=(10^(0.1*Rs)-1)/(10^(0.1*Rp)-1);
-Wpw = tan(pi.*Wp./T); # prewarp
-Wsw = tan(pi.*Ws./T); # prewarp
-
-##pass/stop band to low pass filter transform:
-if (length(Wpw)==2 && length(Wsw)==2)
- wp=1;
- w02 = Wpw(1) * Wpw(2); # Central frequency of stop/pass band (square)
- w3 = w02/Wsw(2);
- w4 = w02/Wsw(1);
- if (w3 > Wsw(1))
- ws = (Wsw(2)-w3)/(Wpw(2)-Wpw(1));
- elseif (w4 < Wsw(2))
- ws = (w4-Wsw(1))/(Wpw(2)-Wpw(1));
- else
- ws = (Wsw(2)-Wsw(1))/(Wpw(2)-Wpw(1));
- endif
-elseif (length(Wpw)==2 && length(Wsw)==1)
- wp=1;
- w02 = Wpw(1) * Wpw(2);
- if (Wsw > Wpw(2))
- w3 = w02/Wsw;
- ws = (Wsw - w3)/(Wpw(2) - Wpw(1));
- else
- w4 = w02/Wsw;
- ws = (w4 - Wsw)/(Wpw(2) - Wpw(1));
- endif
-else
- wp = Wpw;
- ws = Wsw;
-endif
-
-k=wp/ws;
-k1=sqrt(1-k^2);
-q0=(1/2)*((1-sqrt(k1))/(1+sqrt(k1)));
-q= q0 + 2*q0^5 + 15*q0^9 + 150*q0^13; %(....)
-D=(10^(0.1*Rs)-1)/(10^(0.1*Rp)-1);
-
-n=ceil(log10(16*D)/log10(1/q));
+ n=ceil(log10(16*D)/log10(1/q));
+endfunction
Modified: trunk/octave-forge/main/signal/inst/fht.m
===================================================================
--- trunk/octave-forge/main/signal/inst/fht.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/fht.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -12,7 +12,6 @@
##
## You should have received a copy of the GNU General Public License
## along with this program; If not, see <http://www.gnu.org/licenses/>.
-##
## -*- texinfo -*-
## @deftypefn{Function File} {m = } fht ( d, n, dim )
Modified: trunk/octave-forge/main/signal/inst/filtfilt.m
===================================================================
--- trunk/octave-forge/main/signal/inst/filtfilt.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/filtfilt.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -4,7 +4,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -30,18 +30,6 @@
## y = filtfilt(b,a,x); z = filter(b,a,x); % apply filter
## plot(t,x,';data;',t,y,';filtfilt;',t,z,';filter;')
-## Changelog:
-## 2000 02 pki...@us...
-## - pad with zeros to load up the state vector on filter reverse.
-## - add example
-## 2007 12 po...@gn...
-## - use filtic to compute initial and final states
-## - work for multiple columns as well
-## 2008 12 lc...@es...
-## - fixed instability issues with IIR filters and noisy inputs
-## - initial states computed according to Likhterov & Kopeika, 2003
-## - use of a "reflection method" to reduce end effects
-## - added some basic tests
## TODO: (pkienzle) My version seems to have similar quality to matlab,
## but both are pretty bad. They do remove gross lag errors, though.
@@ -49,7 +37,7 @@
function y = filtfilt(b, a, x)
if (nargin != 3)
- usage("y=filtfilt(b,a,x)");
+ print_usage;
end
rotate = (size(x,1)==1);
if rotate, # a row vector
@@ -134,4 +122,3 @@
%! assert (y, [yr.' ys.']);
%! y2 = filtfilt(b.', a.', [r.' s.']);
%! assert (y, y2);
-
Modified: trunk/octave-forge/main/signal/inst/filtic.m
===================================================================
--- trunk/octave-forge/main/signal/inst/filtic.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/filtic.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
@@ -35,12 +35,10 @@
## input vector x and output vector y
##
-## Author: David Billinghurst
-
function zf = filtic(b,a,y,x)
if (nargin>4 || nargin<3) || (nargout>1)
- usage ("zf = filtic (b, a, y [,x])");
+ print_usage;
endif
if nargin < 4, x = []; endif
Modified: trunk/octave-forge/main/signal/inst/fir1.m
===================================================================
--- trunk/octave-forge/main/signal/inst/fir1.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/fir1.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -2,7 +2,7 @@
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
+## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
Modified: trunk/octave-forge/main/signal/inst/fir2.m
===================================================================
--- trunk/octave-forge/main/signal/inst/fir2.m 2012-01-18 08:20:01 UTC (rev 9530)
+++ trunk/octave-forge/main/signal/inst/fir2.m 2012-01-18 14:02:54 UTC (rev 9531)
@@ -1,17 +1,17 @@
## Copyright (C) 2000 Paul Kienzle <pki...@us...>
##
-## This program is free software; you can redistribute it and/or modify
-## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
-## (at your option) any later version.
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
##
-## This program is distributed in the hope that it will be useful,
-## but WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-## GNU General Public License for more details.
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
##
-## You should have received a copy of the GNU General Public License
-## along with this program; If not, see <http://www.gnu.org/licenses/>.
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
## usage: b = fir2(n, f, m [, grid_n [, ramp_n]] [, window])
##
@@ -43,19 +43,10 @@
## [h, w] = freqz(fir2(100,f,m));
## plot(f,m,';target response;',w/pi,abs(h),';filter response;');
-## Feb 27, 2000 PAK
-## use ramping on any transition less than ramp_n units
-## use 2^nextpow2(n+1) for expanded grid size if grid is too small
-## 2001-01-30 PAK
-## set default ramp length to grid_n/20 (i.e., pi/20 radians)
-## use interp1 to interpolate the grid points
-## better(?) handling of 0 and pi frequency points.
-## added some demos
-
function b = fir2(n, f, m, grid_n, ramp_n, window)
if nargin < 3 || nargin > 6
- usage("b = fir2(n, f, m [, grid_n [, ramp_n]] [, window])");
+ print_usage;
endif
## verify frequency and magnitude vectors are reasonable
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2012-03-10 12:40:23
|
Revision: 9805
http://octave.svn.sourceforge.net/octave/?rev=9805&view=rev
Author: jpicarbajal
Date: 2012-03-10 12:40:15 +0000 (Sat, 10 Mar 2012)
Log Message:
-----------
signal: Applying patches from devian packaging group. Generator Rafael Laboissiere ra...@la...
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/filtfilt.m
trunk/octave-forge/main/signal/inst/fir1.m
trunk/octave-forge/main/signal/inst/fir2.m
trunk/octave-forge/main/signal/inst/impinvar.m
trunk/octave-forge/main/signal/inst/invimpinvar.m
trunk/octave-forge/main/signal/inst/residued.m
Modified: trunk/octave-forge/main/signal/inst/filtfilt.m
===================================================================
--- trunk/octave-forge/main/signal/inst/filtfilt.m 2012-03-10 04:18:08 UTC (rev 9804)
+++ trunk/octave-forge/main/signal/inst/filtfilt.m 2012-03-10 12:40:15 UTC (rev 9805)
@@ -105,11 +105,11 @@
%! y = filtfilt(b, a, r+s);
%! assert (size(r), size(y));
%! assert (mean(abs(y)) < 1e3);
-%! assert (corrcoef(s(250:750), y(250:750)) > .95)
+%! assert (corr(s(250:750), y(250:750)) > .95)
%! [b,a] = butter(2, [4e-4 8e-2]);
%! yb = filtfilt(b, a, r+s);
%! assert (mean(abs(yb)) < 1e3);
-%! assert (corrcoef(y, yb) > .99)
+%! assert (corr(y, yb) > .99)
%!test
%! randn('state',0);
Modified: trunk/octave-forge/main/signal/inst/fir1.m
===================================================================
--- trunk/octave-forge/main/signal/inst/fir1.m 2012-03-10 04:18:08 UTC (rev 9804)
+++ trunk/octave-forge/main/signal/inst/fir1.m 2012-03-10 12:40:15 UTC (rev 9805)
@@ -92,7 +92,7 @@
if rem(n,2)==1 && m(2*bands)==1,
warning("n must be even for highpass and bandstop filters. Incrementing.");
n = n+1;
- if isvector(window) && isreal(window)
+ if isvector(window) && isreal(window) && !ischar(window)
## Extend the window using interpolation
M = length(window);
if M == 1,
Modified: trunk/octave-forge/main/signal/inst/fir2.m
===================================================================
--- trunk/octave-forge/main/signal/inst/fir2.m 2012-03-10 04:18:08 UTC (rev 9804)
+++ trunk/octave-forge/main/signal/inst/fir2.m 2012-03-10 12:40:15 UTC (rev 9805)
@@ -72,7 +72,7 @@
if length(ramp_n)>1, w=ramp_n; ramp_n=grid_n/20; endif
if nargin < 6, window=w; endif
if isempty(window), window=hamming(n+1); endif
- if !isreal(window), window=feval(window, n+1); endif
+ if !isreal(window) || ischar(window), window=feval(window, n+1); endif
if length(window) != n+1, usage("window must be of length n+1"); endif
## make sure grid is big enough for the window
Modified: trunk/octave-forge/main/signal/inst/impinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/impinvar.m 2012-03-10 04:18:08 UTC (rev 9804)
+++ trunk/octave-forge/main/signal/inst/impinvar.m 2012-03-10 12:40:15 UTC (rev 9805)
@@ -123,7 +123,7 @@
%!
%! % calculate impulse response of continuous time system
%! % at discrete time intervals 1/fs
-%! ys=impulse(s,1,(n-1)/fs,n);
+%! ys=impulse(s,(n-1)/fs,1/fs)';
%!
%! % impulse response of discrete time system
%! yz=filter(bz,az,[1 zeros(1,n-1)]);
Modified: trunk/octave-forge/main/signal/inst/invimpinvar.m
===================================================================
--- trunk/octave-forge/main/signal/inst/invimpinvar.m 2012-03-10 04:18:08 UTC (rev 9804)
+++ trunk/octave-forge/main/signal/inst/invimpinvar.m 2012-03-10 12:40:15 UTC (rev 9805)
@@ -126,7 +126,7 @@
%!
%! % calculate impulse response of continuous time system
%! % at discrete time intervals 1/fs
-%! ys=impulse(s,1,(n-1)/fs,n);
+%! ys=impulse(s,(n-1)/fs,1/fs)';
%!
%! % impulse response of discrete time system
%! yz=filter(bz,az,[1 zeros(1,n-1)]);
Modified: trunk/octave-forge/main/signal/inst/residued.m
===================================================================
--- trunk/octave-forge/main/signal/inst/residued.m 2012-03-10 04:18:08 UTC (rev 9804)
+++ trunk/octave-forge/main/signal/inst/residued.m 2012-03-10 12:40:15 UTC (rev 9805)
@@ -15,9 +15,9 @@
%% -*- texinfo -*-
%% @deftypefn {Function File} {[@var{r}, @var{p}, @var{f}, @var{m}] =} residued (@var{B}, @var{A})
-%% Compute the partial fraction expansion (PFE) of filter
+%% Compute the partial fraction expansion (PFE) of filter
%% @math{H(z) = B(z)/A(z)}.
-%% In the usual PFE function @code{residuez},
+%% In the usual PFE function @code{residuez},
%% the IIR part (poles @var{p} and residues
%% @var{r}) is driven @emph{in parallel} with the FIR part (@var{f}).
%% In this variant (@code{residued}) the IIR part is driven
@@ -26,7 +26,7 @@
%%
%% INPUTS:
%% @var{B} and @var{A} are vectors specifying the digital filter @math{H(z) = B(z)/A(z)}.
-%% Say @code{help filter} for documentation of the @var{B} and @var{A}
+%% Say @code{help filter} for documentation of the @var{B} and @var{A}
%% filter coefficients.
%%
%% RETURNED:
@@ -42,9 +42,9 @@
%% Say @code{test residued verbose} to see a number of examples.
%% @end example
%%
-%% For the theory of operation, see
+%% For the theory of operation, see
%% @indicateurl{http://ccrma.stanford.edu/~jos/filters/residued.html}
-%%
+%%
%% @seealso{residue residued}
%% @end deftypefn
@@ -60,21 +60,21 @@
%
% This is the same result as returned by RESIDUEZ.
% Otherwise, the FIR part f will be nonempty,
- % and the returned filter is
+ % and the returned filter is
%
% H(z) = f(1) + f(2)/z + f(3)/z^2 + ... + f(nf)/z^M + R(z)/z^M
%
% where R(z) is the parallel one-pole filter bank defined above,
% and M is the order of F(z) = length(f)-1 = nb-na.
- %
+ %
% Note, in particular, that the impulse-response of the parallel
% (complex) one-pole filter bank starts AFTER that of the the FIR part.
% In the result returned by RESIDUEZ, R(z) is not divided by z^M,
% so its impulse response starts at time 0 in parallel with f(n).
%
% J.O. Smith, 9/19/05
-
- if nargin==3,
+
+ if nargin==3,
warning("tolerance ignored");
end
NUM = b(:)';
@@ -91,7 +91,7 @@
if f2, error('f2 not empty as expected'); end
end
-%!test
+%!test
%! B=1; A=[1 -1];
%! [r,p,f,m] = residued(B,A);
%! assert({r,p,f,m},{1,1,[],1},100*eps);
@@ -99,38 +99,38 @@
%! assert({r,p,f,m},{r2,p2,f2,m2},100*eps);
% residuez and residued should be identical when length(B)<length(A)
-%!test
+%!test
%! B=[1 -2 1]; A=[1 -1];
%! [r,p,f,m] = residued(B,A);
%! assert({r,p,f,m},{0,1,[1 -1],1},100*eps);
-%!test
+%!test
%! B=[1 -2 1]; A=[1 -0.5];
%! [r,p,f,m] = residued(B,A);
%! assert({r,p,f,m},{0.25,0.5,[1 -1.5],1},100*eps);
-%!test
+%!test
%! B=1; A=[1 -0.75 0.125];
%! [r,p,f,m] = residued(B,A);
%! [r2,p2,f2,m2] = residuez(B,A);
%! assert({r,p,f,m},{r2,p2,f2,m2},100*eps);
% residuez and residued should be identical when length(B)<length(A)
-%!test
+%!test
%! B=1; A=[1 -2 1];
%! [r,p,f,m] = residued(B,A);
%! [r2,p2,f2,m2] = residuez(B,A);
%! assert({r,p,f,m},{r2,p2,f2,m2},100*eps);
% residuez and residued should be identical when length(B)<length(A)
-%!test
+%!test
%! B=[6,2]; A=[1 -2 1];
%! [r,p,f,m] = residued(B,A);
%! [r2,p2,f2,m2] = residuez(B,A);
%! assert({r,p,f,m},{r2,p2,f2,m2},100*eps);
% residuez and residued should be identical when length(B)<length(A)
-%!test
+%!test
%! B=[1 1 1]; A=[1 -2 1];
%! [r,p,f,m] = residued(B,A);
%! assert(r,[0;3],1e-7);
@@ -138,7 +138,7 @@
%! assert(f,1,100*eps);
%! assert(m,[1;2],100*eps);
-%!test
+%!test
%! B=[2 6 6 2]; A=[1 -2 1];
%! [r,p,f,m] = residued(B,A);
%! assert(r,[8;16],3e-7);
@@ -146,7 +146,7 @@
%! assert(f,[2,10],100*eps);
%! assert(m,[1;2],100*eps);
-%!test
+%!test
%! B=[1,6,2]; A=[1 -2 1];
%! [r,p,f,m] = residued(B,A);
%! assert(r,[-1;9],3e-7);
@@ -154,8 +154,8 @@
%! assert(f,1,100*eps);
%! assert(m,[1;2],100*eps);
-%!test
+%!test
%! B=[1 0 0 0 1]; A=[1 0 0 0 -1];
%! [r,p,f,m] = residued(B,A);
-%! assert({r,p,f,m},{[-1/2;1/2;-j/2;j/2],[-1;1;-j;j],1,[1;1;1;1]},100*eps);
+%! assert({r,p,f,m},{[-j/2;j/2;1/2;-1/2],[-j;j;1;-1],1,[1;1;1;1]},100*eps);
% Verified in maxima: ratsimp(%I/2/(1-%I * d) - %I/2/(1+%I * d)); etc.
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mtm...@us...> - 2012-05-09 03:49:16
|
Revision: 10386
http://octave.svn.sourceforge.net/octave/?rev=10386&view=rev
Author: mtmiller
Date: 2012-05-09 03:49:09 +0000 (Wed, 09 May 2012)
Log Message:
-----------
signal: make window length argument consistent across all functions
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/chebwin.m
trunk/octave-forge/main/signal/inst/flattopwin.m
trunk/octave-forge/main/signal/inst/gausswin.m
trunk/octave-forge/main/signal/inst/kaiser.m
trunk/octave-forge/main/signal/inst/rectwin.m
trunk/octave-forge/main/signal/inst/triang.m
trunk/octave-forge/main/signal/inst/tukeywin.m
Modified: trunk/octave-forge/main/signal/inst/chebwin.m
===================================================================
--- trunk/octave-forge/main/signal/inst/chebwin.m 2012-05-09 01:25:06 UTC (rev 10385)
+++ trunk/octave-forge/main/signal/inst/chebwin.m 2012-05-09 03:49:09 UTC (rev 10386)
@@ -13,9 +13,9 @@
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
-## Usage: chebwin (n, at)
+## Usage: chebwin (L, at)
##
-## Returns the filter coefficients of the n-point Dolph-Chebyshev window
+## Returns the filter coefficients of the L-point Dolph-Chebyshev window
## with at dB of attenuation in the stop-band of the corresponding
## Fourier transform.
##
@@ -31,13 +31,13 @@
##
## The window is described in frequency domain by the expression:
##
-## Cheb(n-1, beta * cos(pi * k/n))
+## Cheb(L-1, beta * cos(pi * k/L))
## W(k) = -------------------------------
-## Cheb(n-1, beta)
+## Cheb(L-1, beta)
##
## with
##
-## beta = cosh(1/(n-1) * acosh(10^(at/20))
+## beta = cosh(1/(L-1) * acosh(10^(at/20))
##
## and Cheb(m,x) denoting the m-th order Chebyshev polynomial calculated
## at the point x.
@@ -48,38 +48,38 @@
##
## See also: kaiser
-function w = chebwin (n, at)
+function w = chebwin (L, at)
if (nargin != 2)
print_usage;
- elseif !(isscalar (n) && (n == round(n)) && (n > 0))
- error ("chebwin: n has to be a positive integer");
+ elseif !(isscalar (L) && (L == round(L)) && (L > 0))
+ error ("chebwin: L has to be a positive integer");
elseif !(isscalar (at) && (at == real (at)))
error ("chebwin: at has to be a real scalar");
endif
- if (n == 1)
+ if (L == 1)
w = 1;
else
# beta calculation
gamma = 10^(-at/20);
- beta = cosh(1/(n-1) * acosh(1/gamma));
+ beta = cosh(1/(L-1) * acosh(1/gamma));
# freq. scale
- k = (0:n-1);
- x = beta*cos(pi*k/n);
+ k = (0:L-1);
+ x = beta*cos(pi*k/L);
# Chebyshev window (freq. domain)
- p = cheb(n-1, x);
+ p = cheb(L-1, x);
# inverse Fourier transform
- if (rem(n,2))
+ if (rem(L,2))
w = real(fft(p));
- M = (n+1)/2;
+ M = (L+1)/2;
w = w(1:M)/w(1);
w = [w(M:-1:2) w]';
else
# half-sample delay (even order)
- p = p.*exp(j*pi/n * (0:n-1));
+ p = p.*exp(j*pi/L * (0:L-1));
w = real(fft(p));
- M = n/2+1;
+ M = L/2+1;
w = w/w(2);
w = [w(M:-1:2) w(2:M)]';
endif
Modified: trunk/octave-forge/main/signal/inst/flattopwin.m
===================================================================
--- trunk/octave-forge/main/signal/inst/flattopwin.m 2012-05-09 01:25:06 UTC (rev 10385)
+++ trunk/octave-forge/main/signal/inst/flattopwin.m 2012-05-09 03:49:09 UTC (rev 10386)
@@ -1,15 +1,15 @@
## Author: Paul Kienzle <pki...@us...> (2004)
## This program is granted to the public domain.
-## flattopwin(n, [periodic|symmetric])
+## flattopwin(L, [periodic|symmetric])
##
## Return the window f(w):
##
## f(w) = 1 - 1.93 cos(2 pi w) + 1.29 cos(4 pi w)
## - 0.388 cos(6 pi w) + 0.0322cos(8 pi w)
##
-## where w = i/(n-1) for i=0:n-1 for a symmetric window, or
-## w = i/n for i=0:n-1 for a periodic window. The default
+## where w = i/(L-1) for i=0:L-1 for a symmetric window, or
+## w = i/L for i=0:L-1 for a periodic window. The default
## is symmetric. The returned window is normalized to a peak
## of 1 at w = 0.5.
##
@@ -23,23 +23,23 @@
## [1] Gade, S; Herlufsen, H; (1987) "Use of weighting functions in DFT/FFT
## analysis (Part I)", Bruel & Kjaer Technical Review No.3.
-function w = flattopwin (n, sym)
+function w = flattopwin (L, sym)
if nargin == 0 || nargin > 2
print_usage;
endif
- divisor = n-1;
+ divisor = L-1;
if nargin > 1
match = strmatch(sym,['periodic';'symmetric']);
if isempty(match),
error("window type must be periodic or symmetric");
elseif match == 1
- divisor = n;
+ divisor = L;
else
- divisor = n-1;
+ divisor = L-1;
endif
endif
- x = 2*pi*[0:n-1]'/divisor;
+ x = 2*pi*[0:L-1]'/divisor;
w = (1-1.93*cos(x)+1.29*cos(2*x)-0.388*cos(3*x)+0.0322*cos(4*x))/4.6402;
endfunction
Modified: trunk/octave-forge/main/signal/inst/gausswin.m
===================================================================
--- trunk/octave-forge/main/signal/inst/gausswin.m 2012-05-09 01:25:06 UTC (rev 10385)
+++ trunk/octave-forge/main/signal/inst/gausswin.m 2012-05-09 03:49:09 UTC (rev 10386)
@@ -13,21 +13,21 @@
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
-## usage: w = gausswin(n, a)
+## usage: w = gausswin(L, a)
##
-## Generate an n-point gaussian window of the given width. Use larger a
-## for a narrow window. Use larger n for a smoother curve.
+## Generate an L-point gaussian window of the given width. Use larger a
+## for a narrow window. Use larger L for a smoother curve.
##
## w = exp ( -(a*x)^2/2 )
##
-## for x = linspace(-(n-1)/n, (n-1)/n, n)
+## for x = linspace(-(L-1)/L, (L-1)/L, L)
-function x = gausswin(n, w)
+function x = gausswin(L, w)
if nargin < 1 || nargin > 2
print_usage;
end
if nargin == 1, w = 2.5; endif
- x = exp ( -0.5 * ( w/n * [ -(n-1) : 2 : n-1 ]' ) .^ 2 );
+ x = exp ( -0.5 * ( w/L * [ -(L-1) : 2 : L-1 ]' ) .^ 2 );
endfunction
Modified: trunk/octave-forge/main/signal/inst/kaiser.m
===================================================================
--- trunk/octave-forge/main/signal/inst/kaiser.m 2012-05-09 01:25:06 UTC (rev 10385)
+++ trunk/octave-forge/main/signal/inst/kaiser.m 2012-05-09 03:49:09 UTC (rev 10386)
@@ -14,36 +14,36 @@
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
-## usage: kaiser (n, beta)
+## usage: kaiser (L, beta)
##
-## Returns the filter coefficients of the n-point Kaiser window with
+## Returns the filter coefficients of the L-point Kaiser window with
## parameter beta.
##
## For the definition of the Kaiser window, see A. V. Oppenheim &
## R. W. Schafer, "Discrete-Time Signal Processing".
##
-## The continuous version of width n centered about x=0 is:
+## The continuous version of width L centered about x=0 is:
##
-## besseli(0, beta * sqrt(1-(2*x/n).^2))
-## k(x) = -------------------------------------, n/2 <= x <= n/2
+## besseli(0, beta * sqrt(1-(2*x/L).^2))
+## k(x) = -------------------------------------, L/2 <= x <= L/2
## besseli(0, beta)
##
## See also: kaiserord
-function w = kaiser (n, beta = 0.5)
+function w = kaiser (L, beta = 0.5)
if (nargin < 1)
print_usage;
- elseif !(isscalar (n) && (n == round (n)) && (n > 0))
- error ("kaiser: n has to be a positive integer");
+ elseif !(isscalar (L) && (L == round (L)) && (L > 0))
+ error ("kaiser: L has to be a positive integer");
elseif !(isscalar (beta) && (beta == real (beta)))
error ("kaiser: beta has to be a real scalar");
endif
- if (n == 1)
+ if (L == 1)
w = 1;
else
- m = n - 1;
+ m = L - 1;
k = (0 : m)';
k = 2 * beta / m * sqrt (k .* (m - k));
w = besseli (0, k) / besseli (0, beta);
Modified: trunk/octave-forge/main/signal/inst/rectwin.m
===================================================================
--- trunk/octave-forge/main/signal/inst/rectwin.m 2012-05-09 01:25:06 UTC (rev 10385)
+++ trunk/octave-forge/main/signal/inst/rectwin.m 2012-05-09 03:49:09 UTC (rev 10386)
@@ -14,12 +14,12 @@
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
-## @deftypefn {Function File} {[@var{w}] =} rectwin(@var{n})
-## Return the filter coefficients of a rectangle window of length N.
+## @deftypefn {Function File} {[@var{w}] =} rectwin(@var{L})
+## Return the filter coefficients of a rectangle window of length L.
## @seealso{hamming, hanning}
## @end deftypefn
-function w = rectwin(n)
+function w = rectwin(L)
if (nargin < 1); print_usage; end
- w = ones(round(n),1);
+ w = ones(round(L),1);
endfunction
Modified: trunk/octave-forge/main/signal/inst/triang.m
===================================================================
--- trunk/octave-forge/main/signal/inst/triang.m 2012-05-09 01:25:06 UTC (rev 10385)
+++ trunk/octave-forge/main/signal/inst/triang.m 2012-05-09 03:49:09 UTC (rev 10386)
@@ -13,20 +13,20 @@
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
-## usage: w = triang (n)
+## usage: w = triang (L)
##
-## Returns the filter coefficients of a triangular window of length n.
+## Returns the filter coefficients of a triangular window of length L.
## Unlike the bartlett window, triang does not go to zero at the edges
-## of the window. For odd n, triang(n) is equal to bartlett(n+2) except
+## of the window. For odd L, triang(L) is equal to bartlett(L+2) except
## for the zeros at the edges of the window.
-function w = triang(n)
+function w = triang(L)
if (nargin != 1)
print_usage;
- elseif (!isscalar(n) || n != fix (n) || n < 1)
- error("triang(n): n has to be an integer > 0");
+ elseif (!isscalar(L) || L != fix (L) || L < 1)
+ error("triang: L has to be an integer > 0");
endif
- w = 1 - abs ([-(n-1):2:(n-1)]' / (n+rem(n,2)));
+ w = 1 - abs ([-(L-1):2:(L-1)]' / (L+rem(L,2)));
endfunction
%!error triang
Modified: trunk/octave-forge/main/signal/inst/tukeywin.m
===================================================================
--- trunk/octave-forge/main/signal/inst/tukeywin.m 2012-05-09 01:25:06 UTC (rev 10385)
+++ trunk/octave-forge/main/signal/inst/tukeywin.m 2012-05-09 03:49:09 UTC (rev 10386)
@@ -14,9 +14,9 @@
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
-## @deftypefn {Function File} {@var{w} =} tukeywin (@var{m}, @var{r})
+## @deftypefn {Function File} {@var{w} =} tukeywin (@var{L}, @var{r})
## Return the filter coefficients of a Tukey window (also known as the
-## cosine-tapered window) of length @var{m}. @var{r} defines the ratio
+## cosine-tapered window) of length @var{L}. @var{r} defines the ratio
## between the constant section and and the cosine section. It has to be
## between 0 and 1. The function returns a Hanning window for @var{r}
## egals 0 and a full box for @var{r} egals 1. By default @var{r} is set
@@ -28,7 +28,7 @@
## Page 67, Equation 38.
## @end deftypefn
-function w = tukeywin (m, r = 1/2)
+function w = tukeywin (L, r = 1/2)
if (nargin < 1 || nargin > 2)
print_usage;
@@ -45,23 +45,23 @@
switch r
case 0,
## full box
- w = ones (m, 1);
+ w = ones (L, 1);
case 1,
## Hanning window
- w = hanning (m);
+ w = hanning (L);
otherwise
## cosine-tapered window
- t = linspace(0,1,m)(1:end/2)';
+ t = linspace(0,1,L)(1:end/2)';
w = (1 + cos(pi*(2*t/r-1)))/2;
- w(floor(r*(m-1)/2)+2:end) = 1;
- w = [w; ones(mod(m,2)); flipud(w)];
+ w(floor(r*(L-1)/2)+2:end) = 1;
+ w = [w; ones(mod(L,2)); flipud(w)];
endswitch
endfunction
%!demo
-%! m = 100;
+%! L = 100;
%! r = 1/3;
-%! w = tukeywin (m, r);
-%! title(sprintf("%d-point Tukey window, R = %d/%d", m, [p, q] = rat(r), q));
+%! w = tukeywin (L, r);
+%! title(sprintf("%d-point Tukey window, R = %d/%d", L, [p, q] = rat(r), q));
%! plot(w);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <car...@us...> - 2012-05-12 13:27:16
|
Revision: 10420
http://octave.svn.sourceforge.net/octave/?rev=10420&view=rev
Author: carandraug
Date: 2012-05-12 13:27:10 +0000 (Sat, 12 May 2012)
Log Message:
-----------
butter pei_tseng_notch: run test silently
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/butter.m
trunk/octave-forge/main/signal/inst/pei_tseng_notch.m
Modified: trunk/octave-forge/main/signal/inst/butter.m
===================================================================
--- trunk/octave-forge/main/signal/inst/butter.m 2012-05-12 13:25:42 UTC (rev 10419)
+++ trunk/octave-forge/main/signal/inst/butter.m 2012-05-12 13:27:10 UTC (rev 10420)
@@ -132,7 +132,7 @@
%! data=[sinetone(5,sf,10,1),sinetone(10,sf,10,1),sinetone(50,sf,10,1),sinetone(200,sf,10,1),sinetone(400,sf,10,1)];
%! [b, a] = butter ( 1, 50 / sf2 );
%! filtered = filter ( b, a, data );
-%! damp_db = 20 * log10 ( max ( filtered ( end - sf : end, : ) ) )
+%! damp_db = 20 * log10 ( max ( filtered ( end - sf : end, : ) ) );
%! assert ( [ damp_db( 4 ) - damp_db( 5 ), damp_db( 1 : 3 ) ], [ 6 0 0 -3 ], off_db )
%!test
@@ -140,7 +140,7 @@
%! data=[sinetone(5,sf,10,1),sinetone(10,sf,10,1),sinetone(50,sf,10,1),sinetone(200,sf,10,1),sinetone(400,sf,10,1)];
%! [b, a] = butter ( 4, 50 / sf2 );
%! filtered = filter ( b, a, data );
-%! damp_db = 20 * log10 ( max ( filtered ( end - sf : end, : ) ) )
+%! damp_db = 20 * log10 ( max ( filtered ( end - sf : end, : ) ) );
%! assert ( [ damp_db( 4 ) - damp_db( 5 ), damp_db( 1 : 3 ) ], [ 24 0 0 -3 ], off_db )
%!test
@@ -148,7 +148,7 @@
%! data=[sinetone(5,sf,10,1),sinetone(10,sf,10,1),sinetone(50,sf,10,1),sinetone(200,sf,10,1),sinetone(400,sf,10,1)];
%! [b, a] = butter ( 1, 50 / sf2, "high" );
%! filtered = filter ( b, a, data );
-%! damp_db = 20 * log10 ( max ( filtered ( end - sf : end, : ) ) )
+%! damp_db = 20 * log10 ( max ( filtered ( end - sf : end, : ) ) );
%! assert ( [ damp_db( 2 ) - damp_db( 1 ), damp_db( 3 : end ) ], [ 6 -3 0 0 ], off_db )
%!test
@@ -156,7 +156,7 @@
%! data=[sinetone(5,sf,10,1),sinetone(10,sf,10,1),sinetone(50,sf,10,1),sinetone(200,sf,10,1),sinetone(400,sf,10,1)];
%! [b, a] = butter ( 4, 50 / sf2, "high" );
%! filtered = filter ( b, a, data );
-%! damp_db = 20 * log10 ( max ( filtered ( end - sf : end, : ) ) )
+%! damp_db = 20 * log10 ( max ( filtered ( end - sf : end, : ) ) );
%! assert ( [ damp_db( 2 ) - damp_db( 1 ), damp_db( 3 : end ) ], [ 24 -3 0 0 ], off_db )
%!demo
Modified: trunk/octave-forge/main/signal/inst/pei_tseng_notch.m
===================================================================
--- trunk/octave-forge/main/signal/inst/pei_tseng_notch.m 2012-05-12 13:25:42 UTC (rev 10419)
+++ trunk/octave-forge/main/signal/inst/pei_tseng_notch.m 2012-05-12 13:27:10 UTC (rev 10420)
@@ -90,7 +90,7 @@
%! data=[sinetone(49,sf,10,1),sinetone(50,sf,10,1),sinetone(51,sf,10,1)];
%! [b, a] = pei_tseng_notch ( 50 / sf2, 2 / sf2 );
%! filtered = filter ( b, a, data );
-%! damp_db = 20 * log10 ( max ( filtered ( end - 1000 : end, : ) ) )
+%! damp_db = 20 * log10 ( max ( filtered ( end - 1000 : end, : ) ) );
%! assert ( damp_db, [ -3 -251.9 -3 ], 0.1 )
%!test
@@ -99,10 +99,9 @@
%! data=[sinetone(49.5,sf,10,1),sinetone(50,sf,10,1),sinetone(50.5,sf,10,1)];
%! [b, a] = pei_tseng_notch ( 50 / sf2, 1 / sf2 );
%! filtered = filter ( b, a, data );
-%! damp_db = 20 * log10 ( max ( filtered ( end - 1000 : end, : ) ) )
+%! damp_db = 20 * log10 ( max ( filtered ( end - 1000 : end, : ) ) );
%! assert ( damp_db, [ -3 -240.4 -3 ], 0.1 )
-
%!demo
%! sf = 800; sf2 = sf/2;
%! data=[[1;zeros(sf-1,1)],sinetone(49,sf,1,1),sinetone(50,sf,1,1),sinetone(51,sf,1,1)];
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2012-06-01 12:43:06
|
Revision: 10550
http://octave.svn.sourceforge.net/octave/?rev=10550&view=rev
Author: jpicarbajal
Date: 2012-06-01 12:42:55 +0000 (Fri, 01 Jun 2012)
Log Message:
-----------
signal: edit to xcorr.m and xcov.m texinfo
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/xcorr.m
trunk/octave-forge/main/signal/inst/xcov.m
Modified: trunk/octave-forge/main/signal/inst/xcorr.m
===================================================================
--- trunk/octave-forge/main/signal/inst/xcorr.m 2012-05-31 13:51:53 UTC (rev 10549)
+++ trunk/octave-forge/main/signal/inst/xcorr.m 2012-06-01 12:42:55 UTC (rev 10550)
@@ -15,60 +15,106 @@
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
-## usage: [R, lag] = xcorr (X [, Y] [, maxlag] [, scale])
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{R}, @var{lag}] =} xcorr ( @var{X} )
+## @deftypefnx {Function File} {@dots{} =} xcorr ( @var{X}, @var{Y} )
+## @deftypefnx {Function File} {@dots{} =} xcorr ( @dots{}, @var{maxlag})
+## @deftypefnx {Function File} {@dots{} =} xcorr ( @dots{}, @var{scale})
+## Estimates the cross-correlation.
##
-## Estimate the cross correlation R_xy(k) of vector arguments X and Y
-## or, if Y is omitted, estimate autocorrelation R_xx(k) of vector X,
-## for a range of lags k specified by argument "maxlag". If X is a
-## matrix, each column of X is correlated with itself and every other
+## Estimate the cross correlation R_xy(k) of vector arguments @var{X} and @var{Y}
+## or, if @var{Y} is omitted, estimate autocorrelation R_xx(k) of vector @var{X},
+## for a range of lags k specified by argument "maxlag". If @var{X} is a
+## matrix, each column of @var{X} is correlated with itself and every other
## column.
##
## The cross-correlation estimate between vectors "x" and "y" (of
-## length N) for lag "k" is given by
-## R_xy(k) = sum_{i=1}^{N}{x_{i+k} conj(y_i),
-## where data not provided (for example x(-1), y(N+1)) is zero.
+## length N) for lag "k" is given by
##
-## ARGUMENTS
-## X [non-empty; real or complex; vector or matrix] data
+## @iftex
+## @tex
+## $$ R_{xy}(k) = \sum_{i=1}^{N} x_{i+k} \conj(y_i),
+## @end tex
+## @end iftex
+## @ifnottex
+## @example
+## N
+## R_xy(k) = sum x_@{i+k@} conj(y_i),
+## i=1
+## @end example
+## @end ifnottex
##
-## Y [real or complex vector] data
-## If X is a matrix (not a vector), Y must be omitted.
-## Y may be omitted if X is a vector; in this case xcorr
-## estimates the autocorrelation of X.
+## where data not provided (for example x(-1), y(N+1)) is zero. Note the
+## definition of cross-correlation given above. To compute a
+## cross-correlation consistent with the field of statistics, see @command{xcov}.
##
-## maxlag [integer scalar] maximum correlation lag
-## If omitted, the default value is N-1, where N is the
-## greater of the lengths of X and Y or, if X is a matrix,
-## the number of rows in X.
+## @strong{ARGUMENTS}
+## @table @var
+## @item X
+## [non-empty; real or complex; vector or matrix] data
##
-## scale [character string] specifies the type of scaling applied
-## to the correlation vector (or matrix). is one of:
-## 'none' return the unscaled correlation, R,
-## 'biased' return the biased average, R/N,
-## 'unbiased' return the unbiassed average, R(k)/(N-|k|),
-## 'coeff' return the correlation coefficient, R/(rms(x).rms(y)),
-## where "k" is the lag, and "N" is the length of X.
-## If omitted, the default value is "none".
-## If Y is supplied but does not have the ame length as X,
-## scale must be "none".
-##
-## RETURNED VARIABLES
-## R array of correlation estimates
-## lag row vector of correlation lags [-maxlag:maxlag]
+## @item Y
+## [real or complex vector] data
##
-## The array of correlation estimates has one of the following forms.
-## (1) Cross-correlation estimate if X and Y are vectors.
-## (2) Autocorrelation estimate if is a vector and Y is omitted,
-## (3) If X is a matrix, R is an matrix containing the cross-
-## correlation estimate of each column with every other column.
-## Lag varies with the first index so that R has 2*maxlag+1
-## rows and P^2 columns where P is the number of columns in X.
-## If Rij(k) is the correlation between columns i and j of X
-## R(k+maxlag+1,P*(i-1)+j) == Rij(k)
-## for lag k in [-maxlag:maxlag], or
-## R(:,P*(i-1)+j) == xcorr(X(:,i),X(:,j)).
-## "reshape(R(k,:),P,P)" is the cross-correlation matrix for X(k,:).
+## If @var{X} is a matrix (not a vector), @var{Y} must be omitted.
+## @var{Y} may be omitted if @var{X} is a vector; in this case xcorr
+## estimates the autocorrelation of @var{X}.
##
+## @item maxlag
+## [integer scalar] maximum correlation lag
+## If omitted, the default value is N-1, where N is the
+## greater of the lengths of @var{X} and @var{Y} or, if @var{X} is a matrix,
+## the number of rows in @var{X}.
+##
+## @item scale
+## [character string] specifies the type of scaling applied
+## to the correlation vector (or matrix). is one of:
+## @table @samp
+## @item none
+## return the unscaled correlation, R,
+## @item biased
+## return the biased average, R/N,
+## @item unbiased
+## return the unbiassed average, R(k)/(N-|k|),
+## @item coeff
+## return the correlation coefficient, R/(rms(x).rms(y)),
+## where "k" is the lag, and "N" is the length of @var{X}.
+## If omitted, the default value is "none".
+## If @var{Y} is supplied but does not have the same length as @var{X},
+## scale must be "none".
+## @end table
+## @end table
+##
+## @strong{RETURNED VARIABLES}
+## @table @var
+## @item R
+## array of correlation estimates
+## @item lag
+## row vector of correlation lags [-maxlag:maxlag]
+## @end table
+##
+## The array of correlation estimates has one of the following forms:
+## (1) Cross-correlation estimate if @var{X} and @var{Y} are vectors.
+##
+## (2) Autocorrelation estimate if is a vector and @var{Y} is omitted.
+##
+## (3) If @var{X} is a matrix, R is an matrix containing the cross-correlation
+## estimate of each column with every other column. Lag varies with the first
+## index so that R has 2*maxlag+1 rows and P^2 columns where P is the number of
+## columns in @var{X}.
+##
+## If Rij(k) is the correlation between columns i and j of @var{X}
+##
+## @code{R(k+maxlag+1,P*(i-1)+j) == Rij(k)}
+##
+## for lag k in [-maxlag:maxlag], or
+##
+## @code{R(:,P*(i-1)+j) == xcorr(X(:,i),X(:,j))}.
+##
+## @code{reshape(R(k,:),P,P)} is the cross-correlation matrix for @code{X(k,:)}.
+##
+## @seealso{xcov}
+## @end deftypefn
## The cross-correlation estimate is calculated by a "spectral" method
## in which the FFT of the first vector is multiplied element-by-element
@@ -89,7 +135,7 @@
## ( hankel(x(1:k),x(k:N-k)) * y ) ./ N
function [R, lags] = xcorr (X, Y, maxlag, scale)
-
+
if (nargin < 1 || nargin > 4)
print_usage;
endif
@@ -110,7 +156,7 @@
endif
## assign defaults to missing arguments
- if isvector(X)
+ if isvector(X)
## if isempty(Y), Y=X; endif ## this line disables code for autocorr'n
N = max(length(X),length(Y));
else
@@ -121,7 +167,7 @@
## check argument values
if isempty(X) || isscalar(X) || ischar(Y) || ! ismatrix(X)
- error("xcorr: X must be a vector or matrix");
+ error("xcorr: X must be a vector or matrix");
endif
if isscalar(Y) || ischar(Y) || (!isempty(Y) && !isvector(Y))
error("xcorr: Y must be a vector");
@@ -130,7 +176,7 @@
error("xcorr: X must be a vector if Y is specified");
endif
if !isscalar(maxlag) || !isreal(maxlag) || maxlag<0 || fix(maxlag)!=maxlag
- error("xcorr: maxlag must be a single non-negative integer");
+ error("xcorr: maxlag must be a single non-negative integer");
endif
##
## sanity check on number of requested lags
@@ -141,7 +187,7 @@
if (maxlag > N-1)
pad_result = maxlag - (N - 1);
maxlag = N - 1;
- %error("xcorr: maxlag must be less than length(X)");
+ %error("xcorr: maxlag must be less than length(X)");
else
pad_result = 0;
endif
@@ -149,15 +195,15 @@
!strcmp(scale,'none')
error("xcorr: scale must be 'none' if length(X) != length(Y)")
endif
-
+
P = columns(X);
M = 2^nextpow2(N + maxlag);
- if !isvector(X)
+ if !isvector(X)
## For matrix X, correlate each column "i" with all other "j" columns
R = zeros(2*maxlag+1,P^2);
## do FFTs of padded column vectors
- pre = fft (postpad (prepad (X, N+maxlag), M) );
+ pre = fft (postpad (prepad (X, N+maxlag), M) );
post = conj (fft (postpad (X, M)));
## do autocorrelations (each column with itself)
@@ -181,7 +227,7 @@
post = fft( postpad(X(:),M) );
cor = ifft( post .* conj(post) );
R = [ conj(cor(maxlag+1:-1:2)) ; cor(1:maxlag+1) ];
- else
+ else
## compute cross-correlation of X and Y
## If one of X & Y is a row vector, the other can be a column vector.
pre = fft( postpad( prepad( X(:), length(X)+maxlag ), M) );
@@ -193,7 +239,7 @@
## if inputs are real, outputs should be real, so ignore the
## insignificant complex portion left over from the FFT
if isreal(X) && (isempty(Y) || isreal(Y))
- R=real(R);
+ R=real(R);
endif
## correct for bias
@@ -218,7 +264,7 @@
elseif !strcmp(scale, 'none')
error("xcorr: scale must be 'biased', 'unbiased', 'coeff' or 'none'");
endif
-
+
## Pad result if necessary
## (most likely is not required, use "if" to avoid uncessary code)
## At this point, lag varies with the first index in R;
@@ -229,7 +275,7 @@
endif
## Correct the shape (transpose) so it is the same as the first input vector
if isvector(X) && P > 1
- R = R.';
+ R = R.';
endif
## return the lag indices if desired
@@ -241,7 +287,7 @@
endfunction
##------------ Use brute force to compute the correlation -------
-##if !isvector(X)
+##if !isvector(X)
## ## For matrix X, compute cross-correlation of all columns
## R = zeros(2*maxlag+1,P^2);
## for i=1:P
@@ -258,7 +304,7 @@
##elseif isempty(Y)
## ## reshape X so that dot product comes out right
## X = reshape(X, 1, N);
-##
+##
## ## compute autocorrelation for 0:maxlag
## R = zeros (2*maxlag + 1, 1);
## for k=0:maxlag
@@ -271,7 +317,7 @@
## ## reshape and pad so X and Y are the same length
## X = reshape(postpad(X,N), 1, N);
## Y = reshape(postpad(Y,N), 1, N)';
-##
+##
## ## compute cross-correlation
## R = zeros (2*maxlag + 1, 1);
## R(maxlag+1) = X*Y;
Modified: trunk/octave-forge/main/signal/inst/xcov.m
===================================================================
--- trunk/octave-forge/main/signal/inst/xcov.m 2012-05-31 13:51:53 UTC (rev 10549)
+++ trunk/octave-forge/main/signal/inst/xcov.m 2012-06-01 12:42:55 UTC (rev 10550)
@@ -13,24 +13,42 @@
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
-## usage: [c, lag] = xcov (X [, Y] [, maxlag] [, scale])
-##
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{R}, @var{lag}] =} xcov ( @var{X} )
+## @deftypefnx {Function File} {@dots{} =} xcov ( @var{X}, @var{Y} )
+## @deftypefnx {Function File} {@dots{} =} xcov ( @dots{}, @var{maxlag})
+## @deftypefnx {Function File} {@dots{} =} xcov ( @dots{}, @var{scale})
## Compute covariance at various lags [=correlation(x-mean(x),y-mean(y))].
##
-## X: input vector
-## Y: if specified, compute cross-covariance between X and Y,
+## @table @var
+## @item X
+## input vector
+## @item Y
+## if specified, compute cross-covariance between X and Y,
## otherwise compute autocovariance of X.
-## maxlag: is specified, use lag range [-maxlag:maxlag],
+## @item maxlag
+## is specified, use lag range [-maxlag:maxlag],
## otherwise use range [-n+1:n-1].
-## Scale:
-## 'biased' for covariance=raw/N,
-## 'unbiased' for covariance=raw/(N-|lag|),
-## 'coeff' for covariance=raw/(covariance at lag 0),
-## 'none' for covariance=raw
-## 'none' is the default.
+## @item scale:
+## @table @samp
+## @item biased
+## for covariance=raw/N,
+## @item unbiased
+## for covariance=raw/(N-|lag|),
+## @item coeff
+## for covariance=raw/(covariance at lag 0),
+## @item none
+## for covariance=raw
+## @item none
+## is the default.
+## @end table
+## @end table
##
-## Returns the covariance for each lag in the range, plus an
+## Returns the covariance for each lag in the range, plus an
## optional vector of lags.
+##
+## @seealso{xcorr}
+## @end deftypefn
function [retval, lags] = xcov (X, Y, maxlag, scale)
@@ -58,5 +76,5 @@
else
[retval, lags] = xcorr(center(X), maxlag, scale);
endif
-
+
endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X