Test notebook for my Mathematica package

Options[$FrontEnd,NotebookBrowseDirectory]

NotebookDirectory[]

copyright Nasser M. Abbasi
Last updated Jan 10, 2017

This notebook tests each function on nma.m, which is package I maintain that contains some useful functions.

SetDirectory[NotebookDirectory[]]
Get["nma.m"]

?nma`*

butterd

?nma`butterd

butterworth digital filter design

test1

Clear[z];
fs=20000 (*hz*);
pass=2000 (* hz*);
stop=3000 (*hz*);
apass=-1 (*db*);
astop=-15 (*db*);
{filterOrder,Ωc,hz,poles,zeros,a,b}=nma`butterd[fs,z,pass,stop,apass,astop,form->"quadratic"];
TraditionalForm@hz

poles

zeros

test2

{filterOrder,Ωc,hz,poles,zeros,a,b}=nma`butterd[fs,z,pass,stop,apass,astop,method->"impulse",form->"quadratic"];
TraditionalForm@hz

poles

zeros

{filterOrder,Ωc,hz,poles,zeros,a,b}=nma`butterd[fs,z,pass,stop,apass,astop,form->"linear"];
TraditionalForm@hz

butterToBandPass

?nma`butterToBandPass

test1

make a normalized butter low pass analog filter of order 4

Clear[s,omega];
{poles,hs}=nma`getButterworthPolynomial[4,1,s];
TraditionalForm@hs

convert to band pass

gamma1=((2 Pi)/(0.25*10^(-3)));
gamma2=((4 Pi)/(0.25*10^(-3)));
hsBandPass=nma`butterToBandPass[hs,gamma1,gamma2,s];
TraditionalForm@hsBandPass

Plot the magnitude spectrums

ftNew=hsBandPass/.s->(I omega);
Plot[Abs[ftNew],{omega,0,2 gamma2},PlotRange->All,AxesOrigin->{0,0},Epilog->{Line[{{gamma1,0},{gamma1,1}}],Line[{{gamma2,0},{gamma2,1}}]}]

butterToBandStop

?nma`butterToBandStop

test1

make a normalized butter low pass analog filter of order 4

Clear[s,omega];
{poles,hs}=nma`getButterworthPolynomial[4,1,s];
TraditionalForm@hs

convert to band stop

gamma1=((2 Pi)/(0.25*10^(-3)));
gamma2=((4 Pi)/(0.25*10^(-3)));
hsBandStop=nma`butterToBandStop[hs,gamma1,gamma2,s];
TraditionalForm@hsBandStop

Plot the magnitude spectrums

ftNew=hsBandStop/.s->(I omega);
Plot[Abs[ftNew],{omega,0,2gamma2},PlotRange->All,AxesOrigin->{0,0},Epilog->{Line[{{gamma1,0},{gamma1,1}}],Line[{{gamma2,0},{gamma2,1}}]}]

butterToHighPass

?nma`butterToHighPass

test1

make a normalized butter low pass analog filter of order 4

Clear[s,omega];
{poles,hs}=nma`getButterworthPolynomial[4,1,s];
TraditionalForm@hs

convert to high pass

gammaCutoff=((2 Pi)/(0.25*10^(-3)))
hsHigh=nma`butterToHighPass[hs,gammaCutoff,s];
TraditionalForm@hsHigh

Plot the magnitude spectrums

ftNew=hsHigh/.s->(I omega);
Plot[Abs[ftNew],{omega,0,4gammaCutoff},PlotRange->All,AxesOrigin->{0,0},Epilog->Line[{{gammaCutoff,0},{gammaCutoff,1}}]]

butterToLowPass

?nma`butterToLowPass

test1

make a normalized butter low pass analog filter of order 4

Clear[s,ω];
{poles,hs}=nma`getButterworthPolynomial[4,1,s];
TraditionalForm@hs

convert the cut off frequency to a new cutoff frequency

Ωc=((2 Pi)/(0.25*10^(-3)));
hsLow=nma`butterToLowPass[hs,Ωc,s];
TraditionalForm@hsLow

Plot the magnitude spectrums

ftNew=hsLow/.s->(I ω);
Plot[Abs[ftNew],{ω,0,5Ωc},PlotRange->All,AxesOrigin->{0,0},Epilog->Line[{{Ωc,0},{Ωc,1}}]]

centerImageData

?nma`centerImageData

img=Import["ExampleData/lena.tif"];
ImageDimensions[img] (*look at image*)

ImageChannels[img]  (*see how many channels*)

data=ImageData[img];(*get data*)
{nRow,nCol,nChannel}=Dimensions[data]

Map[Image[data[[All,All,#]]]&,Range[1,nChannel]](*look at each channel*)

d=data[[All,All,2]];(*get channel 2*)Image[d,ImageSize->400]

(*center Image first*)
d=nma`centerImageData[d];

Dimensions[d]

d[[2,2]]

(*make FFT and look at spectrum and phase*)
fw=Fourier[d,FourierParameters->{1,1}];

fudgeFactor=100;(*adjust for better viewing as needed*)ArrayPlot[(fudgeFactor*Log[1+Abs@fw]),PlotLabel->"Magnitude spectrum",ImageSize->400]
ArrayPlot[Arg[fw],PlotLabel->"Phase spectrum",ImageSize->400]

convertToLCD

?nma`convertToLCD

convertToLCD[{{5/6, 1/4}, {1/7, 1/8}}]

convertToLCD[{{5/6, 1/4}}]

diagonalBlockMatrix

?nma`diagonalBlockMatrix

(block={{1,2},{3,4}})//MatrixForm
nma`diagonalBlockMatrix[block,4];
MatrixForm[Normal[%]]

diagonalDominantMatrixQ

?nma`diagonalDominantMatrixQ

mat={{4,4},{4,4}};
nma`diagonalDominantMatrixQ[mat]

nma`diagonalDominantMatrixQ[mat,strict->True]

displayHurwitz

?nma`displayHurwitz

nma`displayHurwitz[{11 + 9*s + 8*s^2 + 6*s^3 + 3*s^4 + s^5}, s]

fcPulseTrain

?nma`fcPulseTrain

get nth fourier series complex coeff

h=1;width=1;period=2; 
Table[{n,nma`fcPulseTrain[n,h,width,period]},{n,-4,4}]//TableForm

formatBinary

?nma`formatBinary

nma`formatBinary[4,4]

nma`formatBinary[4,8]

nma`formatBinary[127,8]

nma`formatBinary[127,16]

nma`formatBinary[30,16]

generateLTVControllabilityMatrix

?nma`generateLTVControllabilityMatrix

a={{Cos[t],t},{Exp[t],2}};
b={{1},{t}};
r=nma`generateLTVControllabilityMatrix[a,b,t]

Det[r]

%/.t->1

generateLTVObservabilityMatrix

?nma`generateLTVObservabilityMatrix

a={{0,1,t},{t^2,-t,Exp[t]},{1,-2,1}};
c={{1,t,Exp[-t]},{1,1,2}};
(r=nma`generateLTVObservabilityMatrix[a,c,t])//MatrixForm
MatrixRank[r/.t->1]

getButterworthPolynomial

?nma`getButterworthPolynomial

Generate butterworth H(s) given order and cuttoff

test1

Clear[s,form];
order=2;
cuttoff=.2 Pi;
{poles,hs}=nma`getButterworthPolynomial[order,cuttoff,s];
poles
TraditionalForm@hs

test2

{poles,hs}=nma`getButterworthPolynomial[order,cuttoff,s,form->"linear"];
poles
TraditionalForm@hs

test3

{poles,hs}=nma`getButterworthPolynomial[order,cuttoff,s,form->"quadratic"];
poles
TraditionalForm@hs

test4  find normalized Butterworth

{poles,hs}=nma`getButterworthPolynomial[5,1,s,form->"quadratic"];
poles
TraditionalForm@hs

test 5, even

{poles,hs}=nma`getButterworthPolynomial[4,1,s,form->"quadratic"];
poles
TraditionalForm@hs

getControllableForm

?nma`getControllableForm

a={{0,0,1,0},{0,0,0,1},{-2,-1,0,0},{1,-1,0,0}};
MatrixForm[a]

b={{0},{0},{1},{0}};
MatrixForm[b]

{A0,B0}=getControllableForm[a,b];
MatrixForm[A0]

MatrixForm[B0]

getMagnitudeOfPulseTrainFourierCoeff

?nma`getMagnitudeOfPulseTrainFourierCoeff

nma`getMagnitudeOfPulseTrainFourierCoeff[0,1,10,1,5,3]

getMaxPadding

?nma`getMaxPadding

p1=Plot[Sin[x],{x,-Pi,Pi},ImagePadding->30];
p2=Plot[Cos[x],{x,-Pi,Pi},ImagePadding->20];
Grid[{{p2, p1}} /. Graphics[x__] :> 
    Graphics[x, ImagePadding -> getMaxPadding[{p1, p2}]]]

getPhaseOfPulseTrainFourierCoeff

?nma`getPhaseOfPulseTrainFourierCoeff

nma`getPhaseOfPulseTrainFourierCoeff[0,1,10,1,5]

getPositiveAndNegativeTerms

?nma`getPositiveAndNegativeTerms

Clear[c,b,d,h]
expr=3-c+b+d-h
{neg,pos}=nma`getPositiveAndNegativeTerms[expr]

getStateGainVector

?nma`getStateGainVector

a={{0,0,1,0},{0,0,0,1},{-2,1,0,0},{1,-1,0,0}};
b={{0},{0},{1},{0}};
desiredEigenValues={-1+I,-1+I,-1-I,-1-I};
nma`getStateGainVector[a,b,desiredEigenValues,False]

a={{0,0,1,0},{0,0,0,1},{-2,1,0,0},{1,-1,0,0}};
b={{0},{0},{1},{0}};
desiredEigenValues={-1+I,-1+I,-1-I,-1-I};
nma`getStateGainVector[a,b,desiredEigenValues,True]

hurwitz

?nma`hurwitz

nma`hurwitz[11 + 9*s + 8*s^2 + 6*s^3 + 3*s^4 + s^5, s]; 
Grid[%, Frame -> All]

kharitonovPoly

?nma`kharitonovPoly

p=nma`kharitonovPoly[{{11,12},{9,10},{7,8},{5,6},{3,4},{1,2}},s];
TableForm[p]

leadingDet

?nma`leadingDet

nma`hurwitz[11 + 9*s + 8*s^2 + 6*s^3 + 3*s^4 + s^5, s]
nma`leadingDet[%]

linspace

?nma`linspace

nma`linspace[2,3,5]

nma`linspace[2+I,3,5]

nma`linspace[2+I,3,5]//N

nma`linspace[2,3,2]//N

nma`linspace[2,2.1,10]//N

makePulseTrain

?nma`makePulseTrain

delay=.2;period=1;range=5.;tao=.5; h=1;
currentPulses=nma`makePulseTrain[delay,period,range,tao,h];
Framed@Plot[0,{x,-range,range},PlotRange->{{-range,range},{0,1.2}},Epilog->{Thin,Red,currentPulses}]

matrix2DLexOrdering

?nma`matrix2DLexOrdering

(c=Table[0,{i,3},{j,1}])//MatrixForm

{A0,B0}=nma`matrix2DLexOrdering[c];
A0//MatrixForm

B0//MatrixForm

(c=Table[0,{i,3},{j,2}])//MatrixForm
{A0,B0}=nma`matrix2DLexOrdering[c];
A0//MatrixForm

B0//MatrixForm

(c=Table[0,{i,3},{j,3}])//MatrixForm
{A0,B0}=nma`matrix2DLexOrdering[c];
A0//MatrixForm

B0//MatrixForm

matrix2DredBlackPositions

?nma`matrix2DredBlackPositions

nRow=3;
nCol=3;
A0=Table[0,{i,nRow},{j,nCol}];
{r,b}=nma`matrix2DredBlackPositions[nRow,nCol,internal->True];
(A0[[Sequence@@#]]=1)&/@r;
(A0[[Sequence@@#]]=-1)&/@b;
MatrixForm[A0]

{r,b}=nma`matrix2DredBlackPositions[nRow,nCol,internal->False];
(A0[[Sequence@@#]]=1)&/@r;
(A0[[Sequence@@#]]=-1)&/@b;
MatrixForm[A0]

nRow=2;
nCol=5;
A0=Table[0,{i,nRow},{j,nCol}];
{r,b}=nma`matrix2DredBlackPositions[nRow,nCol,internal->False];
(A0[[Sequence@@#]]=1)&/@r;
(A0[[Sequence@@#]]=-1)&/@b;
MatrixForm[A0]

To apply the 2D Laplacian on a matrix using red-black

mask={{0,1,0},{1,0,1},{0,1,0}};
nRow=5; nCol=5;
A0=Table[RandomInteger[10],{i,nRow},{j,nCol}];
MatrixForm[A0]

{r,b}=nma`matrix2DredBlackPositions[nRow,nCol,internal->True];
Table[A0[[Sequence@@r[[i]] ]]=(1./4)Total[mask*  A0[[  r[[i,1]]-1;;r[[i,1]]+1,r[[i,2]]-1;;r[[i,2]]+1  ]],2],{i,1,Length[r]}];
MatrixForm[A0]
Table[A0[[Sequence@@b[[i]] ]]=(1./4)Total[mask*  A0[[  b[[i,1]]-1;;b[[i,1]]+1,b[[i,2]]-1;;b[[i,2]]+1  ]],2],{i,1,Length[b]}];
MatrixForm[A0]

numIt

?nma`numIt

nma`numIt[23.321,5,1,14]
nma`numIt[23.321,5,2,10]

numItse

?nma`numItse

RandomReal[1]
nma`numItse[%,5,3,10]

RandomReal[10]
nma`numItse[%,5,3,10]

plotFourierTransform

?nma`plotFourierTransform

fm=10;
period=1/fm;
x[t_]:=2 Cos[2 Pi fm t]
Plot[x[t],{t,0, 2/10}]
FourierTransform[x[t],t,w,FourierParameters->{0,2Pi}];

Show[First@nma`plotFourierTransform[%,w,0,2 period,-2,2,Large],AxesLabel->{"w (hz)","F(w)"}]

poissonMatrixOnSquareDirichlet

?nma`poissonMatrixOnSquareDirichlet

nma`poissonMatrixOnSquareDirichlet[3,5]//MatrixForm

nma`poissonMatrixOnSquareDirichlet[2,5] //MatrixForm

nma`poissonMatrixOnSquareDirichlet[4,5] //MatrixForm

polarForm

?nma`polarForm

nma`polarForm[3.0 + 2 I]

quantizationLevel

?nma`quantizationLevel

nBits=8;
mp=16;
x=-6.5;

nma`quantizationLevel[-8.7,16,nBits]

snip

?nma`snip

data1 = {{1, 1.1}, {2, 1.5}, {3, 0.9}, {4, 2.3}, {5, 1.1}};
data2 = {{1, 1001.1}, {2, 1001.5}, {3, 1000.9}, {4, 1002.3}, {5, 1001.1}};
p1 = ListPlot[data1, PlotRange -> All, Joined -> True, Mesh -> Full, 
   PlotStyle -> Red, AxesStyle -> {None, snip[1]}, 
   PlotRangePadding -> None, ImagePadding -> 40];

p2 = ListPlot[data2, PlotRange -> All, Joined -> True, Mesh -> Full, 
   PlotStyle -> Blue, Axes -> {False, True}, 
   AxesStyle -> {None, snip[0]}, PlotRangePadding -> None, 
   ImagePadding -> 40];

Column[{p2, p1} /. Graphics[x__]:> Graphics[x, ImagePadding -> getMaxPadding[{p1, p2}],ImageSize -> 400]]

str

?nma`str

test:=Module[{a=4},Print["a=",nma`str[a]]]
test

stripGraphMouseAppearance

?nma`stripGraphMouseAppearance

nma`stripGraphMouseAppearance[
  TreeGraph[{1 -> 2, 1 -> 3, 2 -> 4, 2 -> 5}, 
  VertexLabels -> {5 -> Placed["five", Tooltip]}, ImagePadding -> 20]
 ]

toOffsetBinary

?nma`toOffsetBinary

mp=16;
nBits=8;
x=-8.7
{delta,quantized,code}=nma`toOffsetBinary[x,mp,nBits]
formatBinary[code,nBits]

mp=1;
nBits=4;
x=1/15.
{delta,quantized,code}=nma`toOffsetBinary[x,mp,nBits]
formatBinary[code,nBits]

mp=1;
nBits=4;
x=-100;
{delta,quantized,code}=nma`toOffsetBinary[x,mp,nBits]
formatBinary[code,nBits]

toOnesComplement

?nma`toOnesComplement

mp=10;
nBits=8;
x= N[-6.5]
{delta,quantized,code}=nma`toOnesComplement[x,mp,nBits]
formatBinary[code,nBits]

mp=1;
nBits=4;
x= -1/17.
{delta,quantized,code}=nma`toOnesComplement[x,mp,nBits]
nma`formatBinary[code,nBits]

mp=1;
nBits=4;
x= 13/17.
{delta,quantized,code}=nma`toOnesComplement[x,mp,nBits]
formatBinary[code,nBits]

mp=16;
nBits=8;
x=-8.7;
{delta,quantized,code}=nma`toOnesComplement[x,mp,nBits]
formatBinary[code,nBits]

mp=1;
nBits=4;
x=100;
{delta,quantized,code}=nma`toOnesComplement[x,mp,nBits]
formatBinary[code,nBits]

mp=1;
nBits=4;
x=-100;
{delta,quantized,code}=nma`toOnesComplement[x,mp,nBits]
formatBinary[code,nBits]

delta=1/16.;
data=Table[{n*delta,formatBinary[nma`toOnesComplement[n*delta,1,4][[3]],4]},{n,2^4-1,-2^4,-1}]//TableForm

toSignMagnitude

?nma`toSignMagnitude

mp=16;
nBits=8;
x=-8.7
{delta,quantized,code}=nma`toSignMagnitude[x,mp,nBits]
formatBinary[code,nBits]

mp=10;
nBits=8;
x=-6.5;
{delta,quantized,code}=nma`toSignMagnitude[x,mp,nBits]
formatBinary[code,nBits]

mp=10;
nBits=8;
x=6.5;
{delta,quantized,code}=nma`toSignMagnitude[x,mp,nBits]
formatBinary[code,nBits]

mp=10;
nBits=8;
x=600;
{delta,quantized,code}=nma`toSignMagnitude[x,mp,nBits]
formatBinary[code,nBits]

mp=10;
nBits=8;
x=-600;
{delta,quantized,code}=nma`toSignMagnitude[x,mp,nBits]
formatBinary[code,nBits]

mp=16;
nBits=8;
x=-8.7;
{delta,quantized,code}=nma`toSignMagnitude[x,mp,nBits]
formatBinary[code,nBits]

mp=10;
nBits=8;
x=600;
{delta,quantized,code}=nma`toSignMagnitude[x,mp,nBits]
formatBinary[code,nBits]

toTwoComplement

?nma`toTwoComplement

mp=16;
nBits=8;
x=-8.7
{delta,quantized,code}=nma`toTwoComplement[x,mp,nBits]
formatBinary[code,nBits]

mp=1;
nBits=4;
x=-8.7
{delta,quantized,code}=nma`toTwoComplement[x,mp,nBits]
formatBinary[code,nBits]

Function to obtain Fourier series approximation given fundamental period and complex fourier coefficients 0..N-1

?nma`makePulseTrain

test on a square pulse

h=1;width=1.5;period=4; n=10;delay=0;range=10*period;
Clear[t,x];

currentPulses=nma`makePulseTrain[delay,period,range,width,h];

Show[{Plot[0,{x,0,range},PlotRange->{{0,range},{-(.2 h),1.5 h}},Epilog->{Thin,Red,currentPulses}]}]

Funtion to divide polynomial by powers of variable and keep powers in form (z^(-n))

dividePolyFormattedV1[poly_,var_,order_?(IntegerQ[#]&&Positive[#]&)]:=Module[{},Apart[poly/var^order]];

(*thanks to Carlos Casar de Araujo on mathgroup for the following 2 functions*)

MakeBoxes[Power[x_,y_?(NumberQ[#]&&Negative[#]&)],TraditionalForm]:=SuperscriptBox[MakeBoxes[x,TraditionalForm],MakeBoxes[y,TraditionalForm]];
MakeBoxes[Times[c_,p:Power[x_,y_?(NumberQ[#]&&Negative[#]&)]],TraditionalForm]:=If[OrderedQ[{c,p}],RowBox[{MakeBoxes[c,TraditionalForm],MakeBoxes[p,TraditionalForm]}],RowBox[{MakeBoxes[p,TraditionalForm],MakeBoxes[c,TraditionalForm]}]];

test

Framed@Style[TraditionalForm[dividePolyFormattedV1[1+x,x,3]],14]
TraditionalForm[dividePolyFormattedV1[1+x+x^3,x,3]]
TraditionalForm[dividePolyFormattedV1[1+x+x^2+x^7,x,3]]

Grid[{
{Framed@Style[TraditionalForm[1/z],14]}},Alignment->Center,Spacings->1
]

Function to convert low pass digital filter to high pass using frequency transformation

test

Clear[z,hz,alpha,cutOffFrequency,wp]
hz[z_]:=(0.001836(1+z^-1)^4)/((1-1.5548 z^-1+0.6493 z^-2)(1-1.4996 z^-1 +0.8482 z^-2));

expectedHz[z_]:=(0.02426(1-z^-1)^4)/((1-1.0416 z^-1+0.4019z^-2)(1-0.5561 z^-1+0.7647 z^-2));
cutOffFrequency=.2Pi; (*rad per second*)
wp=.6 Pi; (*rad per second*)
alpha=-((Cos[((wp+cutOffFrequency)/2)])/(Cos[((wp-cutOffFrequency)/2)]));
foundHz[z_]:=hz[(-((z^-1+alpha  )/(1+alpha z^-1)))];
BodePlot[foundHz[z]]

test2

Clear[z,hz,x,expectedHz,foundHz]
hz[z_] :=((0.175 z^-1)/(1-1.318 z^-1+0.493 z^-2));
expectedHz[z_]:=((0.12 z^-2 - 0.258 z^-1 +0.12)/(0.059 z^-2-0.111z^-1+0.327));
cutOffFrequency=.16Pi;
wp=.4Pi;
alpha=-((Cos[((wp+cutOffFrequency)/2)])/(Cos[((wp-cutOffFrequency)/2)]));

foundHz[z_]:=hz[(-((z^-1+alpha  )/(1+alpha z^-1)))];

BodePlot[foundHz[z]]

Function to convert low pass digital filter to bandpass digital filter using frequency transformation

test

Clear[z,hz,alpha,cutOff,w1,w2,k]
hz[z_]:=(0.001836(1+z^-1)^4)/((1-1.5548 z^-1+0.6493 z^-2)(1-1.4996 z^-1 +0.8482 z^-2));

cutOff=.2 Pi; (*rad per second*)
w1=.4Pi;(*rad per second*)
w2=.6 Pi;(*rad per second*)
alpha=((Cos[((w2+w1)/2)])/(Cos[((w2-w1)/2)]));
k=Cot[((w2-w1)/2)]Tan[cutOff/2];
foundHz[z_]:=Module[{},hz[-N[(((z^-2-((2 alpha k)/(k+1))z^-1 +((k-1)/(k+1)) )/(((k-1)/(k+1))z^-2-((2 alpha k)/(k+1))z^-1+1)))]]];


BodePlot[foundHz[z]]

Function to convert low pass digital filter to bandstop digital filter by frequency transformation

test

Clear[z,hz,alpha,cutOff,w1,w2,k]
hz[z_]:=((0.001836(1+z^-1)^4)/((1-1.5548z^-1+0.6493z^-2)(1-1.4996z^-1+0.8482z^-2)));
cutOff=.2Pi;
w1=.4Pi;
w2=.6Pi;
alpha=((Cos[((w2+w1)/2)])/(Cos[((w2-w1)/2)]));
k=Tan[((w2-w1)/2)]Tan[cutOff/2];

foundHz[z_]:=hz[  (z^(-2) - (2*alpha)/(k+1) * z^(-1)+ (1-k)/(k+1))/((1-k)/(k+1)*z^(-2)-2*alpha/(k+1)*z^(-1)+1)]


BodePlot[foundHz[z]]

checkForSingularity

ClearAll[c,a,b,x];
m=nma`checkForSingularity[(c-(a+b+1)x)/(x(1-x)),(-a b)/(x(1-x)),{0,1,Infinity},x];
Grid[Join[{{"point","limit x p(x)","limit x^2 q(x)"}},m],Frame->All]

ClearAll[c,a,b,x];
m=nma`checkForSingularity[-2 x/(1-x^2),(lam+40(1-x^2)-(mu^2)/(1-x^2))/(1-x^2),{-1,1,Infinity},x];
Grid[Join[{{"point","limit x p(x)","limit x^2 q(x)"}},m],Frame->All]