1* Fourier Series Function for SPICE
2* This script is offered here for learning purposes, even if it is outdated
3* and superseeded by the spec function and especially by the much faster fft function.
4* You may use this script in conjunction with e.g. a ringoscillator output (see
5* the ngspice manual, chapter 17).
6
7.control
8  begin
9
10* Variable argc delivers the number of command line parameters given by the user
11* after the 'spectrum' command
12   if ($argc lt 4)
13 	 echo   Error: Too few arguments.
14	 echo   '      'Spectrum produces a plot containing a fourier series transformation of
15	 echo   '      'the specified vectors
16	 echo   usage: spectrum startfreq stop step vec [[vec] ...]
17	 goto bottom
18   end
19
20* Check if vectors 'time' and any input vector(s) are available
21* argv[n] delivers the command line entries after the 'spectrum' command,
22* starting with argv[1]. $argv[4-len] delivers the value of all tokens,
23* starting with postion 4 till the end of the command line
24   if ( time eq time )
25      foreach vec $argv[4-len]
26         if ( $vec eq $vec )
27	     else
28            goto bottom
29	     end
30      end
31   else
32      echo '       'Spectrum can not work without a time vector from a transient analysis.
33      goto bottom
34   end
35
36   * generate a new plot entitled 'scratch', which will hold intermediate
37   * results and will be discarded after their evaluation.
38   set dt=$curplot
39   set title=$curplottitle
40   set curplot=new
41   set scratch=$curplot
42
43   * A vector 'span' is created in the 'scratch' plot to hold the time difference
44   * of the transient simulation. {$dt}.time allows to access the 'time' vector
45   * from the dt plot (which is normally named 'tranx' with x a consecutoive
46   * integer number, depending on the amount of transient simulations already run
47   * in the present job.
48   let span={$dt}.time[length({$dt}.time)-1]-{$dt}.time[0]
49
50* Calculate the number of steps in all of the spectra to be evaluated below
51   if ($argv[3] gt 0.999/span)
52       let fpoints= ( $argv[2] - $argv[1] ) / $argv[3] +1
53	   if (fpoints < 2)
54	       echo frequency start stop or step not correctly specified
55		   goto reset
56	   end
57   else
58       echo   Error: time span is not long enough for a step frequency of $argv[3] Hz
59	   goto reset
60   end
61   let lent = length({$dt}.time)
62   set lent = "$&lent"
63   let nyquist = {$lent}/2/span
64   if ($argv[2] gt nyquist)
65       echo   Error: The nyquist limit is exceeded, try a frequency less than "$&nyquist" Hz
66	   goto reset
67   end
68   set fpoints="$&fpoints"
69
70   * generate a new plot to hold the spectra
71   set curplot=new
72   set spec=$curplot
73   set curplottitle=$title
74   set curplotname='Spectrum Analysis'
75
76* argv[3] is the third agrgument from the input line
77* spectrum 1 1000MEG 10MEG v(out25)
78* that is the delta frequency
79* The fcn vector(n) creates a vector of length n, its elements have
80* the values 0, 1, 2, 3, ..., n-2, n-1. Each element then is multiplied
81* with the frequency step value.
82   let frequency=vector( $fpoints )*$argv[3]
83
84* Add an frequency offset to each element of vector 'frequency'
85* to suppress the (typically) large dc component.
86   dowhile frequency[1] < ( $argv[1] + 1e-9 )
87       let frequency = frequency + $argv[3]
88   end
89
90* For each input vector given on the command line,
91* create a new vector for complex numbers
92   foreach vec $argv[4-len]
93       let $vec = vector( $fpoints ) + j(vector( $fpoints ))
94	   reshape $vec [{$fpoints}]
95   end
96
97* $scratch is a plot for intermediate results, will be destroyed during cleanup
98* $dt is the plot with the original data
99* $spec is a plot for storing the spectrum
100   set curplot=$scratch
101
102   * some test
103   let npers=1
104   let test = span-2/$argv[3] + 1e-9
105   while test > 0
106       let npers = npers + 1
107       let test = test-1/$argv[3]
108   end
109
110   * Do the spectrum calculations
111   let ircle = 2*pi*max(-1,({$dt}.time-{$dt}.time[{$lent}-1])*{$argv[3]}/npers)
112   let win = 1 - cos(ircle)
113   let ircle = npers*ircle
114   let circle = ircle * ({$spec}.frequency[0]/$argv[3] - 1)
115   let k=vector( $fpoints )
116   foreach k $&k
117      let circle = circle + ircle
118      foreach vec $argv[4-len]
119	      let tmp = win*{$dt}.{$vec}
120          let {$spec}.{$vec}[{$k}] = 2*(mean(cos(circle)*tmp),mean(sin(circle)*tmp))
121      end
122   end
123
124* plot (and write) the generated spectrum
125   set curplot = $spec
126   settype frequency frequency
127   foreach vec $argv[4-len]
128      let spectrum = mag({$vec})
129      plot spectrum
130      write specout.out spectrum
131   end
132
133* If you have an oscillator, fimd its frequency
134* as maximum of vector spectrum or goto end (uncomment next line)
135*  goto cleanup
136   set curplot=$scratch
137   let counter = 0
138   let contents = 0
139   let freqmax = 0
140   let spectrum = {$spec}.spectrum
141
142   foreach spectrum $&spectrum
143       if counter > 4
144         if ( contents < $spectrum )
145            let contents = $spectrum
146			set count = "$&counter"
147			let freqmax = {$spec}.frequency[{$count}]
148         end
149      end
150	  let counter = counter + 1
151   end
152
153   echo
154   echo Osc. frequency at "$&freqmax" Hz
155   echo
156   goto cleanup
157
158   label reset
159   set curplot=$dt
160   label cleanup
161   destroy $scratch
162   unset fpoints dt scratch spec vec k title lent
163   label bottom
164
165 end
166