File Coverage

blib/lib/PDL/DSP/Windows.pm
Criterion Covered Total %
statement 399 565 70.6
branch 111 262 42.3
condition 8 30 26.6
subroutine 84 108 77.7
pod 58 97 59.7
total 660 1062 62.1


line stmt bran cond sub pod time code
1             # This file generated by mkwindows.pl
2             package PDL::DSP::Windows;
3             $PDL::DSP::Windows::VERSION = '0.007';
4 2     2   228614 use base 'Exporter';
  2         4  
  2         193  
5 2     2   10 use strict; use warnings;
  2     2   4  
  2         43  
  2         10  
  2         6  
  2         57  
6 2     2   706 use PDL::LiteF;
  2         877  
  2         11  
7 2     2   187755 use PDL::FFT;
  2         16035  
  2         15  
8 2     2   448 use PDL::Math qw( acos cosh acosh );
  2         4  
  2         10  
9 2     2   221 use PDL::Core qw( topdl );
  2         4  
  2         9  
10 2     2   154 use PDL::MatrixOps qw( eigens_sym );
  2         4  
  2         18  
11 2     2   153 use PDL::Options qw( iparse ifhref );
  2         3  
  2         279  
12              
13             eval { require PDL::LinearAlgebra::Special };
14             my $HAVE_LinearAlgebra = 1 if !$@;
15              
16             eval { require PDL::GSLSF::BESSEL; };
17             my $HAVE_BESSEL = 1 if !$@;
18              
19             eval { require PDL::Graphics::Gnuplot; };
20             my $HAVE_GNUPLOT = 1 if !$@;
21              
22             #eval { require PDL::Graphics::PLplot; };
23             #my $HAVE_PLPLOT = 1 if !$@;
24              
25 2     2   1713 use PDL::Constants qw(PI);
  2         23601  
  2         168  
26 2     2   18 use constant TPI => 2 * PI;
  2         2  
  2         25187  
27              
28             our @ISA = qw(Exporter);
29              
30             our @EXPORT_OK = qw( window list_windows chebpoly cos_mult_to_pow cos_pow_to_mult
31             bartlett bartlett_hann bartlett_hann_per bartlett_per blackman blackman_bnh blackman_bnh_per blackman_ex blackman_ex_per blackman_gen blackman_gen3 blackman_gen3_per blackman_gen4 blackman_gen4_per blackman_gen5 blackman_gen5_per blackman_gen_per blackman_harris blackman_harris4 blackman_harris4_per blackman_harris_per blackman_nuttall blackman_nuttall_per blackman_per bohman bohman_per cauchy cauchy_per chebyshev cos_alpha cos_alpha_per cosine cosine_per dpss dpss_per exponential exponential_per flattop flattop_per gaussian gaussian_per hamming hamming_ex hamming_ex_per hamming_gen hamming_gen_per hamming_per hann hann_matlab hann_per hann_poisson hann_poisson_per kaiser kaiser_per lanczos lanczos_per nuttall nuttall1 nuttall1_per nuttall_per parzen parzen_octave parzen_per poisson poisson_per rectangular rectangular_per triangular triangular_per tukey tukey_per welch welch_per);
32              
33             $PDL::onlinedoc->scan(__FILE__) if $PDL::onlinedoc;
34              
35             our %winsubs;
36             our %winpersubs;
37             our %window_definitions;
38              
39             =encoding utf8
40              
41             =head1 NAME
42              
43             PDL::DSP::Windows - Window functions for signal processing
44              
45             =head1 SYNOPSIS
46              
47             use PDL;
48             use PDL::DSP::Windows('window');
49             my $samples = window( 10, 'tukey', { params => .5 });
50              
51             use PDL;
52             use PDL::DSP::Windows;
53             my $win = new PDL::DSP::Windows(10, 'tukey', { params => .5 });
54             print $win->coherent_gain , "\n";
55             $win->plot;
56              
57             =head1 DESCRIPTION
58              
59             This module provides symmetric and periodic (DFT-symmetric)
60             window functions for use in filtering and spectral analysis.
61             It provides a high-level access subroutine
62             L. This functional interface is sufficient for getting the window
63             samples. For analysis and plotting, etc. an object oriented
64             interface is provided. The functional subroutines must be either explicitly exported, or
65             fully qualified. In this document, the word I refers only to the
66             mathematical window functions, while the word I is used to describe
67             code.
68              
69             Window functions are also known as apodization
70             functions or tapering functions. In this module, each of these
71             functions maps a sequence of C<$N> integers to values called
72             a B. (To confuse matters, the word I also has
73             other meanings when describing window functions.)
74             The functions are often named for authors of journal articles.
75             Be aware that across the literature and software,
76             some functions referred to by several different names, and some names
77             refer to several different functions. As a result, the choice
78             of window names is somewhat arbitrary.
79              
80             The L window function requires
81             L. The L window function requires
82             L. But the remaining window functions may
83             be used if these modules are not installed.
84              
85             The most common and easiest usage of this module is indirect, via some
86             higher-level filtering interface, such as L.
87             The next easiest usage is to return a pdl of real-space samples with the subroutine L.
88             Finally, for analyzing window functions, object methods, such as L,
89             L, L are provided.
90              
91             In the following, first the functional interface (non-object oriented) is described in
92             L. Next, the object methods are described in L.
93             Next the low-level subroutines returning samples for each named window
94             are described in L. Finally,
95             some support routines that may be of interest are described in
96             L.
97              
98             =head1 FUNCTIONAL INTERFACE
99              
100             =head2 window
101              
102             $win = window({OPTIONS});
103             $win = window($N,{OPTIONS});
104             $win = window($N,$name,{OPTIONS});
105             $win = window($N,$name,$params,{OPTIONS});
106             $win = window($N,$name,$params,$periodic);
107              
108             Returns an C<$N> point window of type C<$name>.
109             The arguments may be passed positionally in the order
110             C<$N,$name,$params,$periodic>, or they may be passed by
111             name in the hash C.
112              
113             =head3 EXAMPLES
114              
115             # Each of the following return a 100 point symmetric hamming window.
116              
117             $win = window(100);
118             $win = window(100, 'hamming');
119             $win = window(100, { name => 'hamming' );
120             $win = window({ N=> 100, name => 'hamming' );
121              
122             # Each of the following returns a 100 point symmetric hann window.
123              
124             $win = window(100, 'hann');
125             $win = window(100, { name => 'hann' );
126              
127             # Returns a 100 point periodic hann window.
128              
129             $win = window(100, 'hann', { periodic => 1 } );
130              
131             # Returns a 100 point symmetric Kaiser window with alpha=2.
132              
133             $win = window(100, 'kaiser', { params => 2 });
134              
135             =head3 OPTIONS
136              
137             The options follow default PDL::Options rules-- They may be abbreviated,
138             and are case-insensitive.
139              
140             =over
141              
142             =item B
143              
144             (string) name of window function. Default: C.
145             This selects one of the window functions listed below. Note
146             that the suffix '_per', for periodic, may be ommitted. It
147             is specified with the option C<< periodic => 1 >>
148            
149              
150             =item B
151              
152              
153             ref to array of parameter or parameters for the window-function
154             subroutine. Only some window-function subroutines take
155             parameters. If the subroutine takes a single parameter,
156             it may be given either as a number, or a list of one
157             number. For example C<3> or C<[3]>.
158              
159             =item B
160              
161             number of points in window function (the same as the order
162             of the filter) No default value.
163              
164             =item B
165              
166             If value is true, return a periodic rather than a symmetric window function. Default: 0
167             (that is, false. that is, symmetric.)
168              
169             =back
170              
171             =cut
172              
173             sub window {
174 143     143 1 62112 my $win = new PDL::DSP::Windows(@_);
175 143         320 $win->samples();
176             }
177              
178             =head2 list_windows
179              
180             list_windows
181             list_windows STR
182              
183             C prints the names all of the available windows.
184             C prints only the names of windows matching
185             the string C.
186              
187             =cut
188              
189             sub list_windows {
190 0     0 1 0 my ($expr) = @_;
191 0         0 my @match;
192 0 0       0 if ($expr) {
193 0         0 my @alias;
194 0         0 foreach (sort keys %winsubs) {
195 0 0       0 push(@match,$_) , next if /$expr/i;
196 0 0       0 push(@match, $_ . ' (alias ' . $alias[0] . ')') if @alias = grep(/$expr/i,@{$window_definitions{$_}->{alias}});
  0         0  
197             }
198             }
199             else {
200 0         0 @match = sort keys %winsubs;
201             }
202 0         0 print join(', ',@match),"\n";
203             }
204              
205              
206             =head1 METHODS
207              
208             =head2 new
209              
210             =for usage
211              
212             my $win = new PDL::DSP::Windows(ARGS);
213              
214             =for ref
215              
216             Create an instance of a Windows object. If C are given, the instance
217             is initialized. C are interpreted in exactly the
218             same way as arguments the subroutine L.
219              
220             =for example
221              
222             For example:
223              
224             my $win1 = new PDL::DSP::Windows(8,'hann');
225             my $win2 = new PDL::DSP::Windows( { N => 8, name => 'hann' } );
226              
227             =cut
228              
229             sub new {
230 146     146 1 3409 my $proto = shift;
231 146   33     626 my $class = ref($proto) || $proto;
232 146         251 my $self = {};
233 146         268 bless ($self, $class);
234 146 100       539 $self->init(@_) if (@_);
235 146         271 return $self;
236             }
237              
238             =head2 init
239              
240             =for usage
241              
242             $win->init(ARGS);
243              
244             =for ref
245              
246             Initialize (or reinitialize) a Windows object. ARGS are interpreted in exactly the
247             same way as arguments the subroutine L.
248              
249             =for example
250              
251             For example:
252              
253             my $win = new PDL::DSP::Windows(8,'hann');
254             $win->init(10,'hamming');
255              
256             =cut
257              
258             sub init {
259 159     159 1 8897 my $self = shift;
260              
261 159         900 my $opt = new PDL::Options(
262             {
263             name => 'hamming',
264             periodic => 0, # symmetric or periodic
265             N => undef, # order
266             params => undef,
267             });
268 159         7688 my ($N,$name,$params,$periodic);
269 159 100       419 $N = shift unless ref ($_[0]);
270 159 100       389 $name = shift unless ref ($_[0]);
271 159 100       357 $params = shift unless ref ($_[0]) eq 'HASH';
272 159 100       310 $periodic = shift unless ref ($_[0]);
273 159 100       342 my $iopts = @_ ? shift : {};
274 159         470 my $opts = $opt->options($iopts);
275 159 100       22864 $name = $opts->{name} unless $name;
276 159         302 $name =~ s/_per$//;
277 159 100       341 $N = $opts->{N} unless $N;
278 159 100       454 $params = $opts->{params} unless defined $params;
279 159 100 100     503 $params = [$params] if defined $params and not ref $params;
280 159 50       400 $periodic = $opts->{periodic} unless $periodic;
281 159 100       346 my $ws = $periodic ? \%winpersubs : \%winsubs;
282 159 50       489 if ( not exists $ws->{$name}) {
283 0 0       0 my $perstr = $periodic ? 'periodic' : 'symmetric';
284 0         0 barf "window: Unknown $perstr window '$name'.";
285             }
286 159         315 $self->{name} = $name;
287 159         248 $self->{N} = $N;
288 159         218 $self->{periodic} = $periodic;
289 159         234 $self->{params} = $params;
290 159         333 $self->{code} = $ws->{$name};
291 159         280 $self->{samples} = undef;
292 159         665 $self->{modfreqs} = undef;
293 159         973 return $self;
294             }
295              
296             =head2 samples
297              
298             =for usage
299              
300             $win->samples();
301              
302             =for ref
303              
304             Generate and return a reference to the piddle of $N samples for the window C<$win>.
305             This is the real-space representation of the window.
306              
307             The samples are stored in the object C<$win>, but are regenerated
308             every time C is invoked. See the method
309             L below.
310              
311             =for example
312              
313             For example:
314              
315             my $win = new PDL::DSP::Windows(8,'hann');
316             print $win->samples(), "\n";
317              
318             =cut
319              
320             sub samples {
321 159     159 1 220 my $self = shift;
322 159 100       430 my @args = defined $self->{params} ? ($self->{N}, @{$self->{params}} ) : ($self->{N});
  50         129  
323 159         504 $self->{samples} = $self->{code}->(@args);
324             }
325              
326             =head2 modfreqs
327              
328             =for usage
329              
330             $win->modfreqs();
331              
332             =for ref
333              
334             Generate and return a reference to the piddle of the modulus of the
335             fourier transform of the samples for the window C<$win>.
336              
337             These values are stored in the object C<$win>, but are regenerated
338             every time C is invoked. See the method
339             L below.
340              
341             =head3 options
342              
343             =over
344              
345             =item min_bins => MIN
346              
347             This sets the minimum number of frequency bins.
348             Default 1000. If necessary, the piddle of window samples
349             are padded with zeros before the fourier transform is performed.
350              
351             =back
352              
353             =cut
354              
355             sub modfreqs {
356 2     2 1 15 my ($self,$inopts) = @_;
357 2         11 my %opts = iparse( { min_bins => 1000 }, ifhref($inopts));
358 2         376 my $data = $self->get('samples');
359 2         9 my $n = $data->nelem;
360 2 50       10 my $fn = $n > $opts{min_bins} ? 2 * $n : $opts{min_bins};
361 2         4 $n--;
362 2         6 my $freq = zeroes($fn);
363 2         109 $freq->slice("0:$n") .= $data;
364 2         69 PDL::FFT::realfft($freq);
365 2         351 my $real = zeros($freq);
366 2         107 my $img = zeros($freq);
367 2         103 my $mid = ($freq->nelem)/2 - 1;
368 2         4 my $mid1 = $mid + 1;
369 2         10 $real->slice("0:$mid") .= $freq->slice("$mid:0:-1");
370 2         88 $real->slice("$mid1:-1") .= $freq->slice("0:$mid");
371 2         82 $img->slice("0:$mid") .= $freq->slice("-1:$mid1:-1");
372 2         83 $img->slice("$mid1:-1") .= $freq->slice("$mid1:-1");
373 2         152 my $mod = $real**2 + $img**2;
374 2         48 $self->{modfreqs} = $mod;
375 2         21 return $mod;
376             }
377              
378             =head2 get
379              
380             =for usage
381              
382             my $windata = $win->get('samples');
383              
384             =for ref
385              
386             Get an attribute (or list of attributes) of the window C<$win>.
387             If attribute C is requested, then the samples are created with the
388             method L if they don't exist.
389              
390             =for example
391              
392             For example:
393              
394             my $win = new PDL::DSP::Windows(8,'hann');
395             print $win->get('samples'), "\n";
396              
397             =cut
398              
399             sub get {
400 16     16 1 23 my $self = shift;
401 16         20 my @res;
402 16         37 foreach (@_) {
403 16 50 33     110 $self->samples() if $_ eq 'samples' and not defined $self->{samples};
404 16 50 33     36683 $self->freqs() if $_ eq 'modfreqs' and not defined $self->{modfreqs};
405 16         51 push @res, $self->{$_};
406             };
407 16 50       57 return wantarray ? @res : $res[0];
408             }
409              
410             =head2 get_samples
411              
412             =for usage
413              
414             my $windata = $win->get_samples
415              
416             =for ref
417              
418             Return a reference to the pdl of samples for the Window instance C<$win>.
419             The samples will be generated with the method L if and only if
420             they have not yet been generated.
421              
422             =cut
423              
424             sub get_samples {
425 0     0 1 0 my $self = shift;
426 0 0       0 $self->{samples} ? $self->{samples} : $self->samples;
427             }
428              
429             =head2 get_modfreqs
430              
431             =for usage
432              
433             my $winfreqs = $win->get_modfreqs;
434             my $winfreqs = $win->get_modfreqs({OPTS});
435              
436             =for ref
437              
438             Return a reference to the pdl of the frequency response (modulus of the DFT)
439             for the Window instance C<$win>.
440              
441             Options are passed to the method L.
442             The data are created with L
443             if they don't exist. The data are also created even
444             if they already exist if options are supplied. Otherwise
445             the cached data are returned.
446              
447             =head3 options
448              
449             =over
450              
451             =item min_bins => MIN
452              
453             This sets the minimum number of frequency bins. See
454             L. Default 1000.
455              
456             =back
457              
458             =cut
459              
460             sub get_modfreqs {
461 0     0 1 0 my $self = shift;
462 0         0 my ($in_opts) = @_;
463 0 0 0     0 ($self->{modfreqs} and not $in_opts) ? $self->{modfreqs} : $self->modfreqs($in_opts);
464             }
465              
466             =head2 get_params
467              
468             =for usage
469              
470             my $params = $win->get_params
471              
472             =for ref
473              
474             Create a new array containing the parameter values for the instance C<$win>
475             and return a reference to the array.
476             Note that not all window types take parameters.
477              
478             =cut
479              
480             sub get_params {
481 0     0 1 0 my $self = shift;
482 0         0 $self->{params};
483             }
484              
485             sub get_N {
486 0     0 0 0 my $self = shift;
487 0         0 $self->{N};
488             }
489              
490             =head2 get_name
491              
492             =for usage
493              
494             print $win->get_name , "\n";
495              
496             =for ref
497              
498             Return a name suitable for printing associated with the window $win. This is
499             something like the name used in the documentation for the particular
500             window function. This is static data and does not depend on the instance.
501              
502             =cut
503              
504             sub get_name {
505 0     0 1 0 my $self = shift;
506 0         0 my $wd = $window_definitions{$self->{name}};
507 0 0       0 return $wd->{pfn} . ' window' if $wd->{pfn};
508 0 0 0     0 return $wd->{fn} . ' window' if $wd->{fn} and not $wd->{fn} =~ /^\*/;
509 0 0       0 return $wd->{fn} if $wd->{fn};
510 0         0 return ucfirst($self->{name}) . ' window';
511             }
512              
513             sub get_param_names {
514 0     0 0 0 my $self = shift;
515 0         0 my $wd = $window_definitions{$self->{name}};
516 0 0       0 $wd->{params} ? ref($wd->{params}) ? $wd->{params} : [$wd->{params}] : undef;
    0          
517             }
518              
519             sub format_param_vals {
520 0     0 0 0 my $self = shift;
521 0         0 my $p = $self->get('params');
522 0 0       0 return '' unless $p;
523 0         0 my $names = $self->get_param_names;
524 0         0 my @p = @$p;
525 0         0 my @names = @$names;
526 0 0       0 return '' unless $names;
527 0         0 my @s;
528 0         0 map { s/^\$// } @names;
  0         0  
529 0         0 foreach (@p) {
530 0         0 push @s, (shift @names) . ' = ' . $_;
531             }
532 0         0 join(', ', @s);
533             }
534              
535             sub format_plot_param_vals {
536 0     0 0 0 my $self = shift;
537 0         0 my $ps = $self->format_param_vals;
538 0 0       0 return '' unless $ps;
539 0         0 ': ' . $ps;
540             }
541              
542             =head2 plot
543              
544             =for usage
545              
546             $win->plot;
547              
548             =for ref
549              
550             Plot the samples. Currently, only PDL::Graphics::Gnuplot is supported.
551             The default display type is used.
552              
553             =cut
554              
555             sub plot {
556 0     0 1 0 my $self = shift;
557 0 0       0 barf "PDL::DSP::Windows::plot Gnuplot not available!" unless $HAVE_GNUPLOT;
558 0         0 my $w = $self->get('samples');
559 0         0 my $title = $self->get_name() .$self->format_plot_param_vals;
560 0         0 PDL::Graphics::Gnuplot::plot( title => $title, xlabel => 'Time (samples)',
561             ylabel => 'amplitude', $w );
562 0         0 return $self;
563             }
564              
565             =head2 plot_freq
566              
567             =for usage
568              
569             Can be called like this
570              
571             $win->plot_freq;
572              
573              
574             Or this
575              
576             $win->plot_freq( {ordinate => ORDINATE });
577              
578              
579             =for ref
580              
581             Plot the frequency response (magnitude of the DFT of the window samples).
582             The response is plotted in dB, and the frequency
583             (by default) as a fraction of the Nyquist frequency.
584             Currently, only PDL::Graphics::Gnuplot is supported.
585             The default display type is used.
586              
587             =head3 options
588              
589             =over
590              
591             =item coord => COORD
592              
593             This sets the units of frequency of the co-ordinate axis.
594             C must be one of C, for
595             fraction of the nyquist frequency (range C<-1,1>),
596             C, for fraction of the sampling frequncy (range
597             C<-.5,.5>), or C for frequency bin number (range
598             C<0,$N-1>). The default value is C.
599              
600             =item min_bins => MIN
601              
602             This sets the minimum number of frequency bins. See
603             L. Default 1000.
604              
605             =back
606              
607             =cut
608              
609             sub plot_freq {
610 0     0 1 0 my $self = shift;
611 0         0 my $opt = new PDL::Options(
612             {
613             coord => 'nyquist',
614             min_bins => 1000
615             });
616 0 0       0 my $iopts = @_ ? shift : {};
617 0         0 my $opts = $opt->options($iopts);
618 0 0       0 barf "PDL::DSP::Windows::plot Gnuplot not available!" unless $HAVE_GNUPLOT;
619 0         0 my $mf = $self->get_modfreqs({ min_bins => $opts->{min_bins}});
620 0         0 $mf /= $mf->max;
621 0         0 my $param_str = $self->format_plot_param_vals;
622 0         0 my $title = $self->get_name() . $param_str
623             . ', frequency response. ENBW=' . sprintf("%2.3f",$self->enbw);
624 0         0 my $coord = $opts->{coord};
625 0         0 my ($coordinate_range,$xlab);
626 0 0       0 if ($coord eq 'nyquist') {
    0          
    0          
627 0         0 $coordinate_range = 1;
628 0         0 $xlab = 'Fraction of Nyquist frequency';
629             }
630             elsif ($coord eq 'sample') {
631 0         0 $coordinate_range = .5;
632 0         0 $xlab = 'Fraction of sampling freqeuncy';
633             }
634             elsif ($coord eq 'bin') {
635 0         0 $coordinate_range = ($self->get_N)/2;
636 0         0 $xlab = 'bin';
637             }
638             else {
639 0         0 barf "plot_freq: Unknown ordinate unit specification $coord";
640             }
641 0         0 my $coordinates = zeroes($mf)->xlinvals(-$coordinate_range,$coordinate_range);
642 0         0 my $ylab = 'freqeuncy response (dB)';
643 0         0 PDL::Graphics::Gnuplot::plot(title => $title,
644             xmin => -$coordinate_range, xmax => $coordinate_range,
645             xlabel => $xlab, ylabel => $ylab,
646             with => 'line', $coordinates, 20 * log10($mf) );
647 0         0 return $self;
648             }
649              
650             =head2 enbw
651              
652             =for usage
653              
654             $win->enbw;
655              
656             =for ref
657              
658             Compute and return the equivalent noise bandwidth of the window.
659              
660             =cut
661              
662             sub enbw {
663 14     14 1 19 my $self = shift;
664 14         39 my $w = $self->get('samples'); # hmm have to quote samples here
665 14         7878 ($w->nelem) * ($w**2)->sum / ($w->sum)**2;
666             }
667              
668             =head2 coherent_gain
669              
670             =for usage
671              
672             $win->coherent_gain;
673              
674             =for ref
675              
676             Compute and return the coherent gain (the dc gain) of the window.
677             This is just the average of the samples.
678              
679             =cut
680              
681             sub coherent_gain {
682 0     0 1 0 my $self = shift;
683 0         0 my $w = $self->get('samples');
684 0         0 $w->sum / $w->nelem;
685             }
686              
687              
688             =head2 process_gain
689              
690             =for usage
691              
692             $win->coherent_gain;
693              
694             =for ref
695              
696             Compute and return the processing gain (the dc gain) of the window.
697             This is just the multiplicative inverse of the C.
698              
699             =cut
700              
701             sub process_gain {
702 0     0 1 0 my $self = shift;
703 0         0 1/$self->enbw();
704             }
705              
706             # not quite correct for some reason.
707             # Seems like 10*log10(this) / 1.154
708             # gives the correct answer in decibels
709              
710             =head2 scallop_loss
711              
712             =for usage
713              
714             $win->scallop_loss;
715              
716             =for ref
717              
718             **BROKEN**.
719             Compute and return the scalloping loss of the window.
720              
721             =cut
722              
723             sub scallop_loss {
724 0     0 1 0 my ($w) = @_;
725             # my $x = (sequence($w) - ($w->nelem/2)) * (PI/$w->nelem);
726 0         0 my $x = sequence($w) * (PI/$w->nelem);
727 0         0 sqrt( (($w*cos($x))->sum)**2 + (($w*sin($x))->sum)**2 ) /
728             $w->sum;
729             }
730              
731             =head1 WINDOW FUNCTIONS
732              
733             These window-function subroutines return a pdl of $N samples. For most
734             windows, there are a symmetric and a periodic version. The
735             symmetric versions are functions of $N points, uniformly
736             spaced, and taking values from x_lo through x_hi. Here, a
737             periodic function of C< $N > points is equivalent to its
738             symmetric counterpart of C<$N+1> points, with the final
739             point omitted. The name of a periodic window-function subroutine is the
740             same as that for the corresponding symmetric function, except it
741             has the suffix C<_per>. The descriptions below describe the
742             symmetric version of each window.
743              
744             The term 'Blackman-Harris family' is meant to include the Hamming family
745             and the Blackman family. These are functions of sums of cosines.
746              
747             Unless otherwise noted, the arguments in the cosines of all symmetric
748             window functions are multiples of C<$N> numbers uniformly spaced
749             from C<0> through C<2 pi>.
750              
751             =cut
752              
753             sub bartlett {
754 2 50   2 1 8 barf "bartlett: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
755 2         3 my ($N) = @_;
756 2         8 1 - abs (zeroes($N)->xlinvals(-1,1));
757             }
758             $window_definitions{bartlett} = {
759             alias => [ 'fejer'],
760             };
761             $winsubs{bartlett} = \&bartlett;
762              
763             sub bartlett_per {
764 2 50   2 0 17 barf "bartlett: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
765 2         5 my ($N) = @_;
766 2         9 1 - abs (zeroes($N)->xlinvals(-1, (-1+1*($N-1))/$N));
767             }
768             $window_definitions{bartlett} = {
769             alias => [ 'fejer'],
770             };
771             $winpersubs{bartlett}= \&bartlett_per;
772              
773             sub bartlett_hann {
774 3 50   3 1 10 barf "bartlett_hann: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
775 3         6 my ($N) = @_;
776 3         13 0.62 - 0.48 * abs (zeroes($N)->xlinvals(-0.5,0.5)) + 0.38* (cos(zeroes($N)->xlinvals(-(PI),PI)));
777             }
778             $window_definitions{bartlett_hann} = {
779             fn => q!Bartlett-Hann!,
780             alias => [ 'Modified Bartlett-Hann'],
781             };
782             $winsubs{bartlett_hann} = \&bartlett_hann;
783              
784             sub bartlett_hann_per {
785 2 50   2 0 10 barf "bartlett_hann: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
786 2         4 my ($N) = @_;
787 2         8 0.62 - 0.48 * abs (zeroes($N)->xlinvals(-0.5, (-0.5+0.5*($N-1))/$N)) + 0.38* (cos(zeroes($N)->xlinvals(-(PI), (-(PI)+PI*($N-1))/$N)));
788             }
789             $window_definitions{bartlett_hann} = {
790             fn => q!Bartlett-Hann!,
791             alias => [ 'Modified Bartlett-Hann'],
792             };
793             $winpersubs{bartlett_hann}= \&bartlett_hann_per;
794              
795             sub blackman {
796 3 50   3 1 11 barf "blackman: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
797 3         6 my ($N) = @_;
798 3         12 my $cx = (cos(zeroes($N)->xlinvals(0,TPI)));
799              
800 3         2036 (0.34) + ($cx * ((-0.5) + ($cx * (0.16))));
801             }
802             $window_definitions{blackman} = {
803             fn => q!'classic' Blackman!,
804             };
805             $winsubs{blackman} = \&blackman;
806              
807             sub blackman_per {
808 2 50   2 0 8 barf "blackman: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
809 2         5 my ($N) = @_;
810 2         7 my $cx = (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
811              
812 2         333 (0.34) + ($cx * ((-0.5) + ($cx * (0.16))));
813             }
814             $window_definitions{blackman} = {
815             fn => q!'classic' Blackman!,
816             };
817             $winpersubs{blackman}= \&blackman_per;
818              
819             sub blackman_bnh {
820 2 50   2 1 9 barf "blackman_bnh: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
821 2         5 my ($N) = @_;
822 2         7 my $cx = (cos(zeroes($N)->xlinvals(0,TPI)));
823              
824 2         319 (0.3461008) + ($cx * ((-0.4973406) + ($cx * (0.1565586))));
825             }
826             $window_definitions{blackman_bnh} = {
827             pfn => q!Blackman-Harris (bnh)!,
828             fn => q!*An improved version of the 3-term Blackman-Harris window given by Nuttall (Ref 2, p. 89).!,
829             };
830             $winsubs{blackman_bnh} = \&blackman_bnh;
831              
832             sub blackman_bnh_per {
833 2 50   2 0 9 barf "blackman_bnh: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
834 2         6 my ($N) = @_;
835 2         8 my $cx = (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
836              
837 2         423 (0.3461008) + ($cx * ((-0.4973406) + ($cx * (0.1565586))));
838             }
839             $window_definitions{blackman_bnh} = {
840             pfn => q!Blackman-Harris (bnh)!,
841             fn => q!*An improved version of the 3-term Blackman-Harris window given by Nuttall (Ref 2, p. 89).!,
842             };
843             $winpersubs{blackman_bnh}= \&blackman_bnh_per;
844              
845             sub blackman_ex {
846 2 50   2 1 7 barf "blackman_ex: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
847 2         4 my ($N) = @_;
848 2         9 my $cx = (cos(zeroes($N)->xlinvals(0,TPI)));
849              
850 2         318 (0.349742046431642) + ($cx * ((-0.496560619088564) + ($cx * (0.153697334479794))));
851             }
852             $window_definitions{blackman_ex} = {
853             fn => q!'exact' Blackman!,
854             };
855             $winsubs{blackman_ex} = \&blackman_ex;
856              
857             sub blackman_ex_per {
858 2 50   2 0 8 barf "blackman_ex: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
859 2         6 my ($N) = @_;
860 2         7 my $cx = (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
861              
862 2         309 (0.349742046431642) + ($cx * ((-0.496560619088564) + ($cx * (0.153697334479794))));
863             }
864             $window_definitions{blackman_ex} = {
865             fn => q!'exact' Blackman!,
866             };
867             $winpersubs{blackman_ex}= \&blackman_ex_per;
868              
869             sub blackman_gen {
870 2 50   2 1 8 barf "blackman_gen: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
871 2         5 my ($N,$alpha) = @_;
872 2         7 my $cx = (cos(zeroes($N)->xlinvals(0,TPI)));
873              
874 2         317 (.5 - $alpha) + ($cx * ((-.5) + ($cx * ($alpha))));
875             }
876             $window_definitions{blackman_gen} = {
877             pfn => q!General classic Blackman!,
878             fn => q!*A single parameter family of the 3-term Blackman window. !,
879             params => [ '$alpha'],
880             };
881             $winsubs{blackman_gen} = \&blackman_gen;
882              
883             sub blackman_gen_per {
884 2 50   2 0 10 barf "blackman_gen: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
885 2         6 my ($N,$alpha) = @_;
886 2         7 my $cx = (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
887              
888 2         323 (.5 - $alpha) + ($cx * ((-.5) + ($cx * ($alpha))));
889             }
890             $window_definitions{blackman_gen} = {
891             pfn => q!General classic Blackman!,
892             fn => q!*A single parameter family of the 3-term Blackman window. !,
893             params => [ '$alpha'],
894             };
895             $winpersubs{blackman_gen}= \&blackman_gen_per;
896              
897             sub blackman_gen3 {
898 2 50   2 1 8 barf "blackman_gen3: 4 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 4;
899 2         5 my ($N,$a0,$a1,$a2) = @_;
900 2         10 my $cx = (cos(zeroes($N)->xlinvals(0,TPI)));
901              
902 2         363 ($a0 - $a2) + ($cx * ((-$a1) + ($cx * (2*$a2))));
903             }
904             $window_definitions{blackman_gen3} = {
905             fn => q!*The general form of the Blackman family. !,
906             params => [ '$a0','$a1','$a2'],
907             };
908             $winsubs{blackman_gen3} = \&blackman_gen3;
909              
910             sub blackman_gen3_per {
911 2 50   2 0 9 barf "blackman_gen3: 4 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 4;
912 2         5 my ($N,$a0,$a1,$a2) = @_;
913 2         10 my $cx = (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
914              
915 2         375 ($a0 - $a2) + ($cx * ((-$a1) + ($cx * (2*$a2))));
916             }
917             $window_definitions{blackman_gen3} = {
918             fn => q!*The general form of the Blackman family. !,
919             params => [ '$a0','$a1','$a2'],
920             };
921             $winpersubs{blackman_gen3}= \&blackman_gen3_per;
922              
923             sub blackman_gen4 {
924 2 50   2 1 8 barf "blackman_gen4: 5 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 5;
925 2         5 my ($N,$a0,$a1,$a2,$a3) = @_;
926 2         9 my $cx = (cos(zeroes($N)->xlinvals(0,TPI)));
927            
928 2         323 ($a0 - $a2) + ($cx * ((-$a1 + 3 * $a3) + ($cx * (2*$a2 + $cx * (-4*$a3) ))));
929             }
930             $window_definitions{blackman_gen4} = {
931             fn => q!*The general 4-term Blackman-Harris window. !,
932             params => [ '$a0','$a1','$a2','$a3'],
933             };
934             $winsubs{blackman_gen4} = \&blackman_gen4;
935              
936             sub blackman_gen4_per {
937 2 50   2 0 10 barf "blackman_gen4: 5 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 5;
938 2         6 my ($N,$a0,$a1,$a2,$a3) = @_;
939 2         8 my $cx = (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
940            
941 2         334 ($a0 - $a2) + ($cx * ((-$a1 + 3 * $a3) + ($cx * (2*$a2 + $cx * (-4*$a3) ))));
942             }
943             $window_definitions{blackman_gen4} = {
944             fn => q!*The general 4-term Blackman-Harris window. !,
945             params => [ '$a0','$a1','$a2','$a3'],
946             };
947             $winpersubs{blackman_gen4}= \&blackman_gen4_per;
948              
949             sub blackman_gen5 {
950 2 50   2 1 8 barf "blackman_gen5: 6 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 6;
951 2         4 my ($N,$a0,$a1,$a2,$a3,$a4) = @_;
952 2         9 my $cx = (cos(zeroes($N)->xlinvals(0,TPI)));
953            
954 2         325 ($a0 - $a2 + $a4) + ($cx * ((-$a1 + 3 * $a3) + ($cx * (2*$a2 -8*$a4 + $cx * (-4*$a3 +$cx *(8*$a4)) ))));
955             }
956             $window_definitions{blackman_gen5} = {
957             fn => q!*The general 5-term Blackman-Harris window. !,
958             params => [ '$a0','$a1','$a2','$a3','$a4'],
959             };
960             $winsubs{blackman_gen5} = \&blackman_gen5;
961              
962             sub blackman_gen5_per {
963 2 50   2 0 9 barf "blackman_gen5: 6 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 6;
964 2         7 my ($N,$a0,$a1,$a2,$a3,$a4) = @_;
965 2         7 my $cx = (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
966            
967 2         309 ($a0 - $a2 + $a4) + ($cx * ((-$a1 + 3 * $a3) + ($cx * (2*$a2 -8*$a4 + $cx * (-4*$a3 +$cx *(8*$a4)) ))));
968             }
969             $window_definitions{blackman_gen5} = {
970             fn => q!*The general 5-term Blackman-Harris window. !,
971             params => [ '$a0','$a1','$a2','$a3','$a4'],
972             };
973             $winpersubs{blackman_gen5}= \&blackman_gen5_per;
974              
975             sub blackman_harris {
976 2 50   2 1 8 barf "blackman_harris: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
977 2         4 my ($N) = @_;
978 2         7 my $cx = (cos(zeroes($N)->xlinvals(0,TPI)));
979              
980 2         342 (0.343103) + ($cx * ((-0.49755) + ($cx * (0.15844))));
981             }
982             $window_definitions{blackman_harris} = {
983             fn => q!Blackman-Harris!,
984             alias => [ 'Minimum three term (sample) Blackman-Harris'],
985             };
986             $winsubs{blackman_harris} = \&blackman_harris;
987              
988             sub blackman_harris_per {
989 2 50   2 0 9 barf "blackman_harris: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
990 2         5 my ($N) = @_;
991 2         10 my $cx = (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
992              
993 2         363 (0.343103) + ($cx * ((-0.49755) + ($cx * (0.15844))));
994             }
995             $window_definitions{blackman_harris} = {
996             fn => q!Blackman-Harris!,
997             alias => [ 'Minimum three term (sample) Blackman-Harris'],
998             };
999             $winpersubs{blackman_harris}= \&blackman_harris_per;
1000              
1001             sub blackman_harris4 {
1002 4 50   4 1 14 barf "blackman_harris4: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1003 4         8 my ($N) = @_;
1004 4         16 my $cx = (cos(zeroes($N)->xlinvals(0,TPI)));
1005            
1006 4         1969 (0.21747) + ($cx * ((-0.45325) + ($cx * (0.28256 + $cx * (-0.04672) ))));
1007             }
1008             $window_definitions{blackman_harris4} = {
1009             fn => q!minimum (sidelobe) four term Blackman-Harris!,
1010             alias => [ 'Blackman-Harris'],
1011             };
1012             $winsubs{blackman_harris4} = \&blackman_harris4;
1013              
1014             sub blackman_harris4_per {
1015 2 50   2 0 8 barf "blackman_harris4: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1016 2         4 my ($N) = @_;
1017 2         9 my $cx = (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
1018            
1019 2         311 (0.21747) + ($cx * ((-0.45325) + ($cx * (0.28256 + $cx * (-0.04672) ))));
1020             }
1021             $window_definitions{blackman_harris4} = {
1022             fn => q!minimum (sidelobe) four term Blackman-Harris!,
1023             alias => [ 'Blackman-Harris'],
1024             };
1025             $winpersubs{blackman_harris4}= \&blackman_harris4_per;
1026              
1027             sub blackman_nuttall {
1028 3 50   3 1 10 barf "blackman_nuttall: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1029 3         5 my ($N) = @_;
1030 3         12 my $cx = (cos(zeroes($N)->xlinvals(0,TPI)));
1031            
1032 3         475 (0.2269824) + ($cx * ((-0.4572542) + ($cx * (0.273199 + $cx * (-0.0425644) ))));
1033             }
1034             $window_definitions{blackman_nuttall} = {
1035             fn => q!Blackman-Nuttall!,
1036             };
1037             $winsubs{blackman_nuttall} = \&blackman_nuttall;
1038              
1039             sub blackman_nuttall_per {
1040 2 50   2 0 8 barf "blackman_nuttall: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1041 2         5 my ($N) = @_;
1042 2         7 my $cx = (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
1043            
1044 2         309 (0.2269824) + ($cx * ((-0.4572542) + ($cx * (0.273199 + $cx * (-0.0425644) ))));
1045             }
1046             $window_definitions{blackman_nuttall} = {
1047             fn => q!Blackman-Nuttall!,
1048             };
1049             $winpersubs{blackman_nuttall}= \&blackman_nuttall_per;
1050              
1051             sub bohman {
1052 4 50   4 1 14 barf "bohman: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1053 4         8 my ($N) = @_;
1054 4         16 my $x = abs((zeroes($N)->xlinvals(-1,1)));
1055 4         909 (1-$x)*cos(PI*$x) +(1/PI)*sin(PI*$x);
1056             }
1057             $window_definitions{bohman} = {
1058             };
1059             $winsubs{bohman} = \&bohman;
1060              
1061             sub bohman_per {
1062 2 50   2 0 9 barf "bohman: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1063 2         3 my ($N) = @_;
1064 2         7 my $x = abs((zeroes($N)->xlinvals(-1, (-1+1*($N-1))/$N)));
1065 2         289 (1-$x)*cos(PI*$x) +(1/PI)*sin(PI*$x);
1066             }
1067             $window_definitions{bohman} = {
1068             };
1069             $winpersubs{bohman}= \&bohman_per;
1070              
1071             sub cauchy {
1072 3 50   3 1 11 barf "cauchy: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1073 3         7 my ($N,$alpha) = @_;
1074 3         13 1 / (1 + ((zeroes($N)->xlinvals(-1,1)) * $alpha)**2);
1075             }
1076             $window_definitions{cauchy} = {
1077             params => [ '$alpha'],
1078             alias => [ 'Abel','Poisson'],
1079             };
1080             $winsubs{cauchy} = \&cauchy;
1081              
1082             sub cauchy_per {
1083 2 50   2 0 11 barf "cauchy: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1084 2         4 my ($N,$alpha) = @_;
1085 2         10 1 / (1 + ((zeroes($N)->xlinvals(-1, (-1+1*($N-1))/$N)) * $alpha)**2);
1086             }
1087             $window_definitions{cauchy} = {
1088             params => [ '$alpha'],
1089             alias => [ 'Abel','Poisson'],
1090             };
1091             $winpersubs{cauchy}= \&cauchy_per;
1092              
1093             sub chebyshev {
1094 2 50   2 1 6 barf "chebyshev: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1095 2         4 my ($N,$at) = @_;
1096            
1097 2         3 my ($M,$M1,$pos,$pos1);
1098 0         0 my $cw;
1099 2         63 my $beta = cosh(1/($N-1) * acosh(1/(10**(-$at/20))));
1100 2         72 my $k = sequence($N);
1101 2         189 my $x = $beta * cos(PI*$k/$N);
1102 2         114 $cw = chebpoly($N-1,$x);
1103 2 100       8 if ( $N % 2 ) { # odd
1104 1         3 $M1 = ($N+1)/2;
1105 1         3 $M = $M1 - 1;
1106 1         2 $pos = 0;
1107 1         1 $pos1 = 1;
1108 1         5 PDL::FFT::realfft($cw);
1109             }
1110             else { # half-sample delay (even order)
1111 1         4 my $arg = PI/$N * sequence($N);
1112 1         102 my $cw_im = $cw * sin($arg);
1113 1         25 $cw *= cos($arg);
1114 1         27 PDL::FFT::fftnd($cw,$cw_im);
1115 1         113 $M1 = $N/2;
1116 1         2 $M = $M1-1;
1117 1         2 $pos = 1;
1118 1         7 $pos1 = 0;
1119             }
1120 2         133 $cw /= ($cw->at($pos));
1121 2         45 my $cwout = zeroes($N);
1122 2         96 $cwout->slice("0:$M") .= $cw->slice("$M:0:-1");
1123 2         88 $cwout->slice("$M1:-1") .= $cw->slice("$pos1:$M");
1124 2         90 $cwout /= max($cwout);
1125 2         148 $cwout;
1126             ;
1127             }
1128             $window_definitions{chebyshev} = {
1129             params => [ '$at'],
1130             alias => [ 'Dolph-Chebyshev'],
1131             };
1132             $winsubs{chebyshev} = \&chebyshev;
1133              
1134             sub cos_alpha {
1135 5 50   5 1 14 barf "cos_alpha: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1136 5         11 my ($N,$alpha) = @_;
1137 5         17 (sin(zeroes($N)->xlinvals(0,PI)))**$alpha ;
1138             }
1139             $window_definitions{cos_alpha} = {
1140             params => [ '$alpha'],
1141             alias => [ 'Power-of-cosine'],
1142             };
1143             $winsubs{cos_alpha} = \&cos_alpha;
1144              
1145             sub cos_alpha_per {
1146 2 50   2 0 8 barf "cos_alpha: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1147 2         5 my ($N,$alpha) = @_;
1148 2         9 (sin(zeroes($N)->xlinvals(0, PI*($N-1)/$N)))**$alpha ;
1149             }
1150             $window_definitions{cos_alpha} = {
1151             params => [ '$alpha'],
1152             alias => [ 'Power-of-cosine'],
1153             };
1154             $winpersubs{cos_alpha}= \&cos_alpha_per;
1155              
1156             sub cosine {
1157 3 50   3 1 10 barf "cosine: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1158 3         5 my ($N) = @_;
1159 3         11 (sin(zeroes($N)->xlinvals(0,PI)));
1160             }
1161             $window_definitions{cosine} = {
1162             alias => [ 'sine'],
1163             };
1164             $winsubs{cosine} = \&cosine;
1165              
1166             sub cosine_per {
1167 2 50   2 0 7 barf "cosine: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1168 2         4 my ($N) = @_;
1169 2         8 (sin(zeroes($N)->xlinvals(0, PI*($N-1)/$N)));
1170             }
1171             $window_definitions{cosine} = {
1172             alias => [ 'sine'],
1173             };
1174             $winpersubs{cosine}= \&cosine_per;
1175              
1176             sub dpss {
1177 0 0   0 1 0 barf "dpss: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1178 0         0 my ($N,$beta) = @_;
1179            
1180 0 0       0 barf 'dpss: PDL::LinearAlgebra not installed.' unless $HAVE_LinearAlgebra;
1181 0 0 0     0 barf "dpss: $beta not between 0 and $N." unless
1182             $beta >= 0 and $beta <= $N;
1183 0         0 $beta /= ($N/2);
1184 0         0 my $k = sequence($N);
1185 0         0 my $s = sin(PI*$beta*$k)/$k;
1186 0         0 $s->slice('0') .= $beta;
1187 0         0 my ($ev,$e) = eigens_sym(PDL::LinearAlgebra::Special::mtoeplitz($s));
1188 0         0 my $i = $e->maximum_ind;
1189 0         0 $ev->slice("($i)")->copy;
1190             ;
1191             }
1192             $window_definitions{dpss} = {
1193             fn => q!Digital Prolate Spheroidal Sequence (DPSS)!,
1194             params => [ '$beta'],
1195             alias => [ 'sleppian'],
1196             };
1197             $winsubs{dpss} = \&dpss;
1198              
1199             sub dpss_per {
1200 0 0   0 0 0 barf "dpss: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1201 0         0 my ($N,$beta) = @_;
1202 0         0 $N++;
1203            
1204 0 0       0 barf 'dpss: PDL::LinearAlgebra not installed.' unless $HAVE_LinearAlgebra;
1205 0 0 0     0 barf "dpss: $beta not between 0 and $N." unless
1206             $beta >= 0 and $beta <= $N;
1207 0         0 $beta /= ($N/2);
1208 0         0 my $k = sequence($N);
1209 0         0 my $s = sin(PI*$beta*$k)/$k;
1210 0         0 $s->slice('0') .= $beta;
1211 0         0 my ($ev,$e) = eigens_sym(PDL::LinearAlgebra::Special::mtoeplitz($s));
1212 0         0 my $i = $e->maximum_ind;
1213 0         0 $ev->slice("($i),0:-2")->copy;
1214             ;
1215             }
1216             $window_definitions{dpss} = {
1217             fn => q!Digital Prolate Spheroidal Sequence (DPSS)!,
1218             params => [ '$beta'],
1219             alias => [ 'sleppian'],
1220             };
1221             $winpersubs{dpss}= \&dpss_per;
1222              
1223             sub exponential {
1224 2 50   2 1 8 barf "exponential: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1225 2         5 my ($N) = @_;
1226 2         7 2 ** (1 - abs (zeroes($N)->xlinvals(-1,1))) - 1;
1227             }
1228             $window_definitions{exponential} = {
1229             };
1230             $winsubs{exponential} = \&exponential;
1231              
1232             sub exponential_per {
1233 2 50   2 0 9 barf "exponential: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1234 2         3 my ($N) = @_;
1235 2         9 2 ** (1 - abs (zeroes($N)->xlinvals(-1, (-1+1*($N-1))/$N))) - 1;
1236             }
1237             $window_definitions{exponential} = {
1238             };
1239             $winpersubs{exponential}= \&exponential_per;
1240              
1241             sub flattop {
1242 4 50   4 1 13 barf "flattop: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1243 4         7 my ($N) = @_;
1244 4         16 my $cx = (cos(zeroes($N)->xlinvals(0,TPI)));
1245            
1246 4         1951 (-0.05473684) + ($cx * ((-0.165894739) + ($cx * (0.498947372 + $cx * (-0.334315788 +$cx *(0.055578944)) ))));
1247             }
1248             $window_definitions{flattop} = {
1249             fn => q!flat top!,
1250             };
1251             $winsubs{flattop} = \&flattop;
1252              
1253             sub flattop_per {
1254 2 50   2 0 8 barf "flattop: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1255 2         4 my ($N) = @_;
1256 2         7 my $cx = (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
1257            
1258 2         317 (-0.05473684) + ($cx * ((-0.165894739) + ($cx * (0.498947372 + $cx * (-0.334315788 +$cx *(0.055578944)) ))));
1259             }
1260             $window_definitions{flattop} = {
1261             fn => q!flat top!,
1262             };
1263             $winpersubs{flattop}= \&flattop_per;
1264              
1265             sub gaussian {
1266 2 50   2 1 10 barf "gaussian: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1267 2         3 my ($N,$beta) = @_;
1268 2         9 exp (-0.5 * ($beta * (zeroes($N)->xlinvals(-1,1)) )**2);
1269             }
1270             $window_definitions{gaussian} = {
1271             params => [ '$beta'],
1272             alias => [ 'Weierstrass'],
1273             };
1274             $winsubs{gaussian} = \&gaussian;
1275              
1276             sub gaussian_per {
1277 2 50   2 0 9 barf "gaussian: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1278 2         3 my ($N,$beta) = @_;
1279 2         9 exp (-0.5 * ($beta * (zeroes($N)->xlinvals(-1, (-1+1*($N-1))/$N)) )**2);
1280             }
1281             $window_definitions{gaussian} = {
1282             params => [ '$beta'],
1283             alias => [ 'Weierstrass'],
1284             };
1285             $winpersubs{gaussian}= \&gaussian_per;
1286              
1287             sub hamming {
1288 6 50   6 1 17 barf "hamming: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1289 6         10 my ($N) = @_;
1290 6         21 0.54 + -0.46 * (cos(zeroes($N)->xlinvals(0,TPI)));
1291             }
1292             $window_definitions{hamming} = {
1293             };
1294             $winsubs{hamming} = \&hamming;
1295              
1296             sub hamming_per {
1297 2 50   2 0 7 barf "hamming: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1298 2         4 my ($N) = @_;
1299 2         9 0.54 + -0.46 * (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
1300             }
1301             $window_definitions{hamming} = {
1302             };
1303             $winpersubs{hamming}= \&hamming_per;
1304              
1305             sub hamming_ex {
1306 2 50   2 1 10 barf "hamming_ex: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1307 2         5 my ($N) = @_;
1308 2         36 0.53836 + -0.46164 * (cos(zeroes($N)->xlinvals(0,TPI)));
1309             }
1310             $window_definitions{hamming_ex} = {
1311             fn => q!'exact' Hamming!,
1312             };
1313             $winsubs{hamming_ex} = \&hamming_ex;
1314              
1315             sub hamming_ex_per {
1316 2 50   2 0 6 barf "hamming_ex: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1317 2         5 my ($N) = @_;
1318 2         9 0.53836 + -0.46164 * (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
1319             }
1320             $window_definitions{hamming_ex} = {
1321             fn => q!'exact' Hamming!,
1322             };
1323             $winpersubs{hamming_ex}= \&hamming_ex_per;
1324              
1325             sub hamming_gen {
1326 2 50   2 1 8 barf "hamming_gen: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1327 2         5 my ($N,$a) = @_;
1328 2         9 $a + -(1-$a) * (cos(zeroes($N)->xlinvals(0,TPI)));
1329             }
1330             $window_definitions{hamming_gen} = {
1331             fn => q!general Hamming!,
1332             params => [ '$a'],
1333             };
1334             $winsubs{hamming_gen} = \&hamming_gen;
1335              
1336             sub hamming_gen_per {
1337 2 50   2 0 8 barf "hamming_gen: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1338 2         4 my ($N,$a) = @_;
1339 2         10 $a + -(1-$a) * (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
1340             }
1341             $window_definitions{hamming_gen} = {
1342             fn => q!general Hamming!,
1343             params => [ '$a'],
1344             };
1345             $winpersubs{hamming_gen}= \&hamming_gen_per;
1346              
1347             sub hann {
1348 5 50   5 1 18 barf "hann: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1349 5         10 my ($N) = @_;
1350 5         21 0.5 + -0.5 * (cos(zeroes($N)->xlinvals(0,TPI)));
1351             }
1352             $window_definitions{hann} = {
1353             alias => [ 'hanning'],
1354             };
1355             $winsubs{hann} = \&hann;
1356              
1357             sub hann_per {
1358 2 50   2 0 10 barf "hann: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1359 2         4 my ($N) = @_;
1360 2         9 0.5 + -0.5 * (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
1361             }
1362             $window_definitions{hann} = {
1363             alias => [ 'hanning'],
1364             };
1365             $winpersubs{hann}= \&hann_per;
1366              
1367             sub hann_matlab {
1368 1 50   1 1 6 barf "hann_matlab: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1369 1         2 my ($N) = @_;
1370 1         5 0.5 - 0.5 * (cos(zeroes($N)->xlinvals(TPI/($N+1),TPI *$N /($N+1))));
1371             }
1372             $window_definitions{hann_matlab} = {
1373             pfn => q!Hann (matlab)!,
1374             fn => q!*Equivalent to the Hann window of N+2 points, with the endpoints (which are both zero) removed.!,
1375             };
1376             $winsubs{hann_matlab} = \&hann_matlab;
1377              
1378             sub hann_poisson {
1379 1 50   1 1 6 barf "hann_poisson: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1380 1         2 my ($N,$alpha) = @_;
1381 1         5 0.5 * (1 + (cos(zeroes($N)->xlinvals(-(PI),PI)))) * exp (-$alpha * abs (zeroes($N)->xlinvals(-1,1)));
1382             }
1383             $window_definitions{hann_poisson} = {
1384             fn => q!Hann-Poisson!,
1385             params => [ '$alpha'],
1386             };
1387             $winsubs{hann_poisson} = \&hann_poisson;
1388              
1389             sub hann_poisson_per {
1390 0 0   0 0 0 barf "hann_poisson: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1391 0         0 my ($N,$alpha) = @_;
1392 0         0 0.5 * (1 + (cos(zeroes($N)->xlinvals(-(PI), (-(PI)+PI*($N-1))/$N)))) * exp (-$alpha * abs (zeroes($N)->xlinvals(-1, (-1+1*($N-1))/$N)));
1393             }
1394             $window_definitions{hann_poisson} = {
1395             fn => q!Hann-Poisson!,
1396             params => [ '$alpha'],
1397             };
1398             $winpersubs{hann_poisson}= \&hann_poisson_per;
1399              
1400             sub kaiser {
1401 0 0   0 1 0 barf "kaiser: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1402 0         0 my ($N,$beta) = @_;
1403            
1404 0 0       0 barf "kaiser: PDL::GSLSF not installed" unless $HAVE_BESSEL;
1405 0         0 $beta *= PI;
1406 0         0 my @n = PDL::GSLSF::BESSEL::gsl_sf_bessel_In ($beta * sqrt(1 - (zeroes($N)->xlinvals(-1,1)) **2),0);
1407 0         0 my @d = PDL::GSLSF::BESSEL::gsl_sf_bessel_In($beta,0);
1408 0         0 (shift @n)/(shift @d);
1409             }
1410             $window_definitions{kaiser} = {
1411             params => [ '$beta'],
1412             alias => [ 'Kaiser-Bessel'],
1413             };
1414             $winsubs{kaiser} = \&kaiser;
1415              
1416             sub kaiser_per {
1417 0 0   0 0 0 barf "kaiser: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1418 0         0 my ($N,$beta) = @_;
1419            
1420 0 0       0 barf "kaiser: PDL::GSLSF not installed" unless $HAVE_BESSEL;
1421 0         0 $beta *= PI;
1422 0         0 my @n = PDL::GSLSF::BESSEL::gsl_sf_bessel_In ($beta * sqrt(1 - (zeroes($N)->xlinvals(-1, (-1+1*($N-1))/$N)) **2),0);
1423 0         0 my @d = PDL::GSLSF::BESSEL::gsl_sf_bessel_In($beta,0);
1424 0         0 (shift @n)/(shift @d);
1425             }
1426             $window_definitions{kaiser} = {
1427             params => [ '$beta'],
1428             alias => [ 'Kaiser-Bessel'],
1429             };
1430             $winpersubs{kaiser}= \&kaiser_per;
1431              
1432             sub lanczos {
1433 3 50   3 1 11 barf "lanczos: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1434 3         9 my ($N) = @_;
1435            
1436 3         11 my $x = PI * (zeroes($N)->xlinvals(-1,1));
1437 3         807 my $res = sin($x)/$x;
1438 3         1020 my $mid;
1439 3 100       20 $mid = int($N/2), $res->slice($mid) .= 1 if $N % 2;
1440 3         61 $res;;
1441             }
1442             $window_definitions{lanczos} = {
1443             alias => [ 'sinc'],
1444             };
1445             $winsubs{lanczos} = \&lanczos;
1446              
1447             sub lanczos_per {
1448 2 50   2 0 9 barf "lanczos: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1449 2         3 my ($N) = @_;
1450            
1451 2         9 my $x = PI * (zeroes($N)->xlinvals(-1, (-1+1*($N-1))/$N));
1452 2         363 my $res = sin($x)/$x;
1453 2         64 my $mid;
1454 2 100       11 $mid = int($N/2), $res->slice($mid) .= 1 unless $N % 2;
1455 2         49 $res;;
1456             }
1457             $window_definitions{lanczos} = {
1458             alias => [ 'sinc'],
1459             };
1460             $winpersubs{lanczos}= \&lanczos_per;
1461              
1462             sub nuttall {
1463 2 50   2 1 8 barf "nuttall: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1464 2         5 my ($N) = @_;
1465 2         7 my $cx = (cos(zeroes($N)->xlinvals(0,TPI)));
1466            
1467 2         332 (0.2269824) + ($cx * ((-0.4572542) + ($cx * (0.273199 + $cx * (-0.0425644) ))));
1468             }
1469             $window_definitions{nuttall} = {
1470             };
1471             $winsubs{nuttall} = \&nuttall;
1472              
1473             sub nuttall_per {
1474 2 50   2 0 7 barf "nuttall: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1475 2         4 my ($N) = @_;
1476 2         8 my $cx = (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
1477            
1478 2         301 (0.2269824) + ($cx * ((-0.4572542) + ($cx * (0.273199 + $cx * (-0.0425644) ))));
1479             }
1480             $window_definitions{nuttall} = {
1481             };
1482             $winpersubs{nuttall}= \&nuttall_per;
1483              
1484             sub nuttall1 {
1485 0 0   0 1 0 barf "nuttall1: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1486 0         0 my ($N) = @_;
1487 0         0 my $cx = (cos(zeroes($N)->xlinvals(0,TPI)));
1488            
1489 0         0 (0.211536) + ($cx * ((-0.449584) + ($cx * (0.288464 + $cx * (-0.050416) ))));
1490             }
1491             $window_definitions{nuttall1} = {
1492             pfn => q!Nuttall (v1)!,
1493             fn => q!*A window referred to as the Nuttall window.!,
1494             };
1495             $winsubs{nuttall1} = \&nuttall1;
1496              
1497             sub nuttall1_per {
1498 0 0   0 0 0 barf "nuttall1: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1499 0         0 my ($N) = @_;
1500 0         0 my $cx = (cos(zeroes($N)->xlinvals(0, TPI*($N-1)/$N)));
1501            
1502 0         0 (0.211536) + ($cx * ((-0.449584) + ($cx * (0.288464 + $cx * (-0.050416) ))));
1503             }
1504             $window_definitions{nuttall1} = {
1505             pfn => q!Nuttall (v1)!,
1506             fn => q!*A window referred to as the Nuttall window.!,
1507             };
1508             $winpersubs{nuttall1}= \&nuttall1_per;
1509              
1510             sub parzen {
1511 3 50   3 1 11 barf "parzen: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1512 3         5 my ($N) = @_;
1513            
1514 3         43 my $x = zeroes($N)->xlinvals(-1,1);
1515 3         842 my $x1 = $x->where($x <= -.5);
1516 3         572 my $x2 = $x->where( ($x < .5) & ($x > -.5) );
1517 3         305 my $x3 = $x->where($x >= .5);
1518 3         230 $x1 .= 2 * (1-abs($x1))**3;
1519             # $x3 .= 2 * (1-abs($x3))**3;
1520 3         1063 $x3 .= $x1->slice('-1:0:-1');
1521 3         410 $x2 .= 1 - 6 * $x2**2 *(1-abs($x2));
1522 3         480 return $x;
1523             }
1524             $window_definitions{parzen} = {
1525             alias => [ 'Jackson','Valle-Poussin'],
1526             };
1527             $winsubs{parzen} = \&parzen;
1528              
1529             sub parzen_per {
1530 2 50   2 0 8 barf "parzen: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1531 2         4 my ($N) = @_;
1532            
1533 2         8 my $x = zeroes($N)->xlinvals(-1,(-1 + ($N-1))/($N));
1534 2         346 my $x1 = $x->where($x <= -.5);
1535 2         146 my $x2 = $x->where( ($x < .5) & ($x > -.5) );
1536 2         120 my $x3 = $x->where($x >= .5);
1537 2         88 $x1 .= 2 * (1-abs($x1))**3;
1538             # $x3 .= 2 * (1-abs($x3))**3;
1539 2         150 $x3 .= $x1->slice('-1:1:-1');
1540 2         91 $x2 .= 1 - 6 * $x2**2 *(1-abs($x2));
1541 2         183 return $x;
1542             }
1543             $window_definitions{parzen} = {
1544             alias => [ 'Jackson','Valle-Poussin'],
1545             };
1546             $winpersubs{parzen}= \&parzen_per;
1547              
1548             sub parzen_octave {
1549 0 0   0 1 0 barf "parzen_octave: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1550 0         0 my ($N) = @_;
1551            
1552 0         0 my $L = $N-1;
1553 0         0 my $r = ($L/2);
1554 0         0 my $r4 = ($r/2);
1555 0         0 my $n = sequence(2*$r+1)-$r;
1556 0         0 my $n1 = $n->where(abs($n) <= $r4);
1557 0         0 my $n2 = $n->where($n > $r4);
1558 0         0 my $n3 = $n->where($n < -$r4);
1559 0         0 $n1 .= 1 -6.*(abs($n1)/($N/2))**2 + 6*(abs($n1)/($N/2))**3;
1560 0         0 $n2 .= 2.*(1-abs($n2)/($N/2))**3;
1561 0         0 $n3 .= 2.*(1-abs($n3)/($N/2))**3;
1562 0         0 $n;
1563             ;
1564             }
1565             $window_definitions{parzen_octave} = {
1566             fn => q!Parzen!,
1567             };
1568             $winsubs{parzen_octave} = \&parzen_octave;
1569              
1570             sub poisson {
1571 3 50   3 1 10 barf "poisson: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1572 3         6 my ($N,$alpha) = @_;
1573 3         12 exp (-$alpha * abs (zeroes($N)->xlinvals(-1,1)));
1574             }
1575             $window_definitions{poisson} = {
1576             params => [ '$alpha'],
1577             };
1578             $winsubs{poisson} = \&poisson;
1579              
1580             sub poisson_per {
1581 2 50   2 0 7 barf "poisson: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1582 2         5 my ($N,$alpha) = @_;
1583 2         8 exp (-$alpha * abs (zeroes($N)->xlinvals(-1, (-1+1*($N-1))/$N)));
1584             }
1585             $window_definitions{poisson} = {
1586             params => [ '$alpha'],
1587             };
1588             $winpersubs{poisson}= \&poisson_per;
1589              
1590             sub rectangular {
1591 4 50   4 1 16 barf "rectangular: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1592 4         9 my ($N) = @_;
1593 4         48 ones($N);
1594             }
1595             $window_definitions{rectangular} = {
1596             alias => [ 'dirichlet','boxcar'],
1597             };
1598             $winsubs{rectangular} = \&rectangular;
1599              
1600             sub rectangular_per {
1601 2 50   2 0 8 barf "rectangular: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1602 2         4 my ($N) = @_;
1603 2         8 ones($N);
1604             }
1605             $window_definitions{rectangular} = {
1606             alias => [ 'dirichlet','boxcar'],
1607             };
1608             $winpersubs{rectangular}= \&rectangular_per;
1609              
1610             sub triangular {
1611 4 50   4 1 16 barf "triangular: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1612 4         7 my ($N) = @_;
1613 4         18 1 - abs (zeroes($N)->xlinvals(-($N-1)/$N,($N-1)/$N));
1614             }
1615             $window_definitions{triangular} = {
1616             };
1617             $winsubs{triangular} = \&triangular;
1618              
1619             sub triangular_per {
1620 2 50   2 0 9 barf "triangular: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1621 2         4 my ($N) = @_;
1622 2         9 1 - abs (zeroes($N)->xlinvals(-$N/($N+1),-1/($N+1)+($N-1)/($N+1)));
1623             }
1624             $window_definitions{triangular} = {
1625             };
1626             $winpersubs{triangular}= \&triangular_per;
1627              
1628             sub tukey {
1629 4 50   4 1 13 barf "tukey: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1630 4         8 my ($N,$alpha) = @_;
1631            
1632 4 50 33     26 barf("tukey: alpha must be between 0 and 1") unless
1633             $alpha >=0 and $alpha <= 1;
1634 4 50       12 return ones($N) if $alpha == 0;
1635 4         17 my $x = zeroes($N)->xlinvals(0,1);
1636 4         941 my $x1 = $x->where($x < $alpha/2);
1637 4         668 my $x2 = $x->where( ($x <= 1-$alpha/2) & ($x >= $alpha/2) );
1638 4         377 my $x3 = $x->where($x > 1 - $alpha/2);
1639 4         298 $x1 .= 0.5 * ( 1 + cos( PI * (2*$x1/$alpha -1)));
1640 4         692 $x2 .= 1;
1641 4         238 $x3 .= $x1->slice('-1:0:-1');
1642 4         250 return $x;
1643             }
1644             $window_definitions{tukey} = {
1645             params => [ '$alpha'],
1646             alias => [ 'tapered cosine'],
1647             };
1648             $winsubs{tukey} = \&tukey;
1649              
1650             sub tukey_per {
1651 2 50   2 0 8 barf "tukey: 2 arguments expected. Got " . scalar(@_) . ' arguments.' unless @_ == 2;
1652 2         3 my ($N,$alpha) = @_;
1653            
1654 2 50 33     14 barf("tukey: alpha must be between 0 and 1") unless
1655             $alpha >=0 and $alpha <= 1;
1656 2 50       8 return ones($N) if $alpha == 0;
1657 2         8 my $x = zeroes($N)->xlinvals(0,($N-1)/$N);
1658 2         323 my $x1 = $x->where($x < $alpha/2);
1659 2         167 my $x2 = $x->where( ($x <= 1-$alpha/2) & ($x >= $alpha/2) );
1660 2         131 my $x3 = $x->where($x > 1 - $alpha/2);
1661 2         98 $x1 .= 0.5 * ( 1 + cos( PI * (2*$x1/$alpha -1)));
1662 2         230 $x2 .= 1;
1663 2         47 $x3 .= $x1->slice('-1:1:-1');
1664 2         92 return $x;
1665             }
1666             $window_definitions{tukey} = {
1667             params => [ '$alpha'],
1668             alias => [ 'tapered cosine'],
1669             };
1670             $winpersubs{tukey}= \&tukey_per;
1671              
1672             sub welch {
1673 3 50   3 1 12 barf "welch: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1674 3         5 my ($N) = @_;
1675 3         13 1 - (zeroes($N)->xlinvals(-1,1))**2;
1676             }
1677             $window_definitions{welch} = {
1678             alias => [ 'Riez','Bochner','Parzen','parabolic'],
1679             };
1680             $winsubs{welch} = \&welch;
1681              
1682             sub welch_per {
1683 2 50   2 0 10 barf "welch: 1 argument expected. Got " . scalar(@_) . ' arguments.' unless @_ == 1;
1684 2         6 my ($N) = @_;
1685 2         9 1 - (zeroes($N)->xlinvals(-1, (-1+1*($N-1))/$N))**2;
1686             }
1687             $window_definitions{welch} = {
1688             alias => [ 'Riez','Bochner','Parzen','parabolic'],
1689             };
1690             $winpersubs{welch}= \&welch_per;
1691              
1692              
1693             =head1 Symmetric window functions
1694              
1695              
1696              
1697             =head2 bartlett($N)
1698              
1699             The Bartlett window. (Ref 1). Another name for this window is the fejer window. This window is defined by
1700              
1701             1 - abs arr,
1702              
1703             where the points in arr range from -1 through 1.
1704             See also L.
1705              
1706              
1707             =cut
1708            
1709              
1710             =head2 bartlett_hann($N)
1711              
1712             The Bartlett-Hann window. Another name for this window is the Modified Bartlett-Hann window. This window is defined by
1713              
1714             0.62 - 0.48 * abs arr + 0.38* arr1,
1715              
1716             where the points in arr range from -1/2 through 1/2, and arr1 are the cos of points ranging from -PI through PI.
1717              
1718              
1719             =cut
1720            
1721              
1722             =head2 blackman($N)
1723              
1724             The 'classic' Blackman window. (Ref 1). One of the Blackman-Harris family, with coefficients
1725              
1726             a0 = 0.42, a1 = 0.5, a2 = 0.08
1727              
1728              
1729              
1730             =cut
1731            
1732              
1733             =head2 blackman_bnh($N)
1734              
1735             The Blackman-Harris (bnh) window. An improved version of the 3-term Blackman-Harris window given by Nuttall (Ref 2, p. 89). One of the Blackman-Harris family, with coefficients
1736              
1737             a0 = 0.4243801, a1 = 0.4973406, a2 = 0.0782793
1738              
1739              
1740              
1741             =cut
1742            
1743              
1744             =head2 blackman_ex($N)
1745              
1746             The 'exact' Blackman window. (Ref 1). One of the Blackman-Harris family, with coefficients
1747              
1748             a0 = 0.426590713671539, a1 = 0.496560619088564, a2 = 0.0768486672398968
1749              
1750              
1751              
1752             =cut
1753            
1754              
1755             =head2 blackman_gen($N,$alpha)
1756              
1757             The General classic Blackman window. A single parameter family of the 3-term Blackman window. This window is defined by
1758              
1759             my $cx = arr;
1760              
1761             (.5 - $alpha) + ($cx * ((-.5) + ($cx * ($alpha)))),
1762              
1763             where the points in arr are the cos of points ranging from 0 through 2PI.
1764              
1765              
1766             =cut
1767            
1768              
1769             =head2 blackman_gen3($N,$a0,$a1,$a2)
1770              
1771             The general form of the Blackman family. One of the Blackman-Harris family, with coefficients
1772              
1773             a0 = $a0, a1 = $a1, a2 = $a2
1774              
1775              
1776              
1777             =cut
1778            
1779              
1780             =head2 blackman_gen4($N,$a0,$a1,$a2,$a3)
1781              
1782             The general 4-term Blackman-Harris window. One of the Blackman-Harris family, with coefficients
1783              
1784             a0 = $a0, a1 = $a1, a2 = $a2, a3 = $a3
1785              
1786              
1787              
1788             =cut
1789            
1790              
1791             =head2 blackman_gen5($N,$a0,$a1,$a2,$a3,$a4)
1792              
1793             The general 5-term Blackman-Harris window. One of the Blackman-Harris family, with coefficients
1794              
1795             a0 = $a0, a1 = $a1, a2 = $a2, a3 = $a3, a4 = $a4
1796              
1797              
1798              
1799             =cut
1800            
1801              
1802             =head2 blackman_harris($N)
1803              
1804             The Blackman-Harris window. (Ref 1). One of the Blackman-Harris family, with coefficients
1805              
1806             a0 = 0.422323, a1 = 0.49755, a2 = 0.07922
1807              
1808             Another name for this window is the Minimum three term (sample) Blackman-Harris window.
1809              
1810             =cut
1811            
1812              
1813             =head2 blackman_harris4($N)
1814              
1815             The minimum (sidelobe) four term Blackman-Harris window. (Ref 1). One of the Blackman-Harris family, with coefficients
1816              
1817             a0 = 0.35875, a1 = 0.48829, a2 = 0.14128, a3 = 0.01168
1818              
1819             Another name for this window is the Blackman-Harris window.
1820              
1821             =cut
1822            
1823              
1824             =head2 blackman_nuttall($N)
1825              
1826             The Blackman-Nuttall window. One of the Blackman-Harris family, with coefficients
1827              
1828             a0 = 0.3635819, a1 = 0.4891775, a2 = 0.1365995, a3 = 0.0106411
1829              
1830              
1831              
1832             =cut
1833            
1834              
1835             =head2 bohman($N)
1836              
1837             The Bohman window. (Ref 1). This window is defined by
1838              
1839             my $x = abs(arr);
1840             (1-$x)*cos(PI*$x) +(1/PI)*sin(PI*$x),
1841              
1842             where the points in arr range from -1 through 1.
1843              
1844              
1845             =cut
1846            
1847              
1848             =head2 cauchy($N,$alpha)
1849              
1850             The Cauchy window. (Ref 1). Other names for this window are: Abel, Poisson. This window is defined by
1851              
1852             1 / (1 + (arr * $alpha)**2),
1853              
1854             where the points in arr range from -1 through 1.
1855              
1856              
1857             =cut
1858            
1859              
1860             =head2 chebyshev($N,$at)
1861              
1862             The Chebyshev window. The frequency response of this window has C<$at> dB of attenuation in the stop-band.
1863             Another name for this window is the Dolph-Chebyshev window. No periodic version of this window is defined.
1864             This routine gives the same result as the routine B in Octave 3.6.2.
1865              
1866              
1867             =cut
1868            
1869              
1870             =head2 cos_alpha($N,$alpha)
1871              
1872             The Cos_alpha window. (Ref 1). Another name for this window is the Power-of-cosine window. This window is defined by
1873              
1874             arr**$alpha ,
1875              
1876             where the points in arr are the sin of points ranging from 0 through PI.
1877              
1878              
1879             =cut
1880            
1881              
1882             =head2 cosine($N)
1883              
1884             The Cosine window. Another name for this window is the sine window. This window is defined by
1885              
1886             arr,
1887              
1888             where the points in arr are the sin of points ranging from 0 through PI.
1889              
1890              
1891             =cut
1892            
1893              
1894             =head2 dpss($N,$beta)
1895              
1896             The Digital Prolate Spheroidal Sequence (DPSS) window. The parameter C<$beta> is the half-width of the mainlobe, measured in frequency bins. This window maximizes the power in the mainlobe for given C<$N> and C<$beta>.
1897             Another name for this window is the sleppian window.
1898              
1899             =cut
1900            
1901              
1902             =head2 exponential($N)
1903              
1904             The Exponential window. This window is defined by
1905              
1906             2 ** (1 - abs arr) - 1,
1907              
1908             where the points in arr range from -1 through 1.
1909              
1910              
1911             =cut
1912            
1913              
1914             =head2 flattop($N)
1915              
1916             The flat top window. One of the Blackman-Harris family, with coefficients
1917              
1918             a0 = 0.21557895, a1 = 0.41663158, a2 = 0.277263158, a3 = 0.083578947, a4 = 0.006947368
1919              
1920              
1921              
1922             =cut
1923            
1924              
1925             =head2 gaussian($N,$beta)
1926              
1927             The Gaussian window. (Ref 1). Another name for this window is the Weierstrass window. This window is defined by
1928              
1929             exp (-0.5 * ($beta * arr )**2),
1930              
1931             where the points in arr range from -1 through 1.
1932              
1933              
1934             =cut
1935            
1936              
1937             =head2 hamming($N)
1938              
1939             The Hamming window. (Ref 1). One of the Blackman-Harris family, with coefficients
1940              
1941             a0 = 0.54, a1 = 0.46
1942              
1943              
1944              
1945             =cut
1946            
1947              
1948             =head2 hamming_ex($N)
1949              
1950             The 'exact' Hamming window. (Ref 1). One of the Blackman-Harris family, with coefficients
1951              
1952             a0 = 0.53836, a1 = 0.46164
1953              
1954              
1955              
1956             =cut
1957            
1958              
1959             =head2 hamming_gen($N,$a)
1960              
1961             The general Hamming window. (Ref 1). One of the Blackman-Harris family, with coefficients
1962              
1963             a0 = $a, a1 = (1-$a)
1964              
1965              
1966              
1967             =cut
1968            
1969              
1970             =head2 hann($N)
1971              
1972             The Hann window. (Ref 1). One of the Blackman-Harris family, with coefficients
1973              
1974             a0 = 0.5, a1 = 0.5
1975              
1976             Another name for this window is the hanning window. See also L.
1977              
1978              
1979             =cut
1980            
1981              
1982             =head2 hann_matlab($N)
1983              
1984             The Hann (matlab) window. Equivalent to the Hann window of N+2 points, with the endpoints (which are both zero) removed. No periodic version of this window is defined.
1985             This window is defined by
1986              
1987             0.5 - 0.5 * arr,
1988              
1989             where the points in arr are the cosine of points ranging from 2PI/($N+1) through 2PI*$N/($N+1).
1990             This routine gives the same result as the routine B in Matlab.
1991             See also L.
1992              
1993              
1994             =cut
1995            
1996              
1997             =head2 hann_poisson($N,$alpha)
1998              
1999             The Hann-Poisson window. (Ref 1). This window is defined by
2000              
2001             0.5 * (1 + arr1) * exp (-$alpha * abs arr),
2002              
2003             where the points in arr range from -1 through 1, and arr1 are the cos of points ranging from -PI through PI.
2004              
2005              
2006             =cut
2007            
2008              
2009             =head2 kaiser($N,$beta)
2010              
2011             The Kaiser window. (Ref 1). The parameter C<$beta> is the approximate half-width of the mainlobe, measured in frequency bins.
2012             Another name for this window is the Kaiser-Bessel window. This window is defined by
2013              
2014            
2015             barf "kaiser: PDL::GSLSF not installed" unless $HAVE_BESSEL;
2016             $beta *= PI;
2017             my @n = PDL::GSLSF::BESSEL::gsl_sf_bessel_In ($beta * sqrt(1 - arr **2),0);
2018             my @d = PDL::GSLSF::BESSEL::gsl_sf_bessel_In($beta,0);
2019             (shift @n)/(shift @d),
2020              
2021             where the points in arr range from -1 through 1.
2022              
2023              
2024             =cut
2025            
2026              
2027             =head2 lanczos($N)
2028              
2029             The Lanczos window. Another name for this window is the sinc window. This window is defined by
2030              
2031            
2032             my $x = PI * arr;
2033             my $res = sin($x)/$x;
2034             my $mid;
2035             $mid = int($N/2), $res->slice($mid) .= 1 if $N % 2;
2036             $res;,
2037              
2038             where the points in arr range from -1 through 1.
2039              
2040              
2041             =cut
2042            
2043              
2044             =head2 nuttall($N)
2045              
2046             The Nuttall window. One of the Blackman-Harris family, with coefficients
2047              
2048             a0 = 0.3635819, a1 = 0.4891775, a2 = 0.1365995, a3 = 0.0106411
2049              
2050             See also L.
2051              
2052              
2053             =cut
2054            
2055              
2056             =head2 nuttall1($N)
2057              
2058             The Nuttall (v1) window. A window referred to as the Nuttall window. One of the Blackman-Harris family, with coefficients
2059              
2060             a0 = 0.355768, a1 = 0.487396, a2 = 0.144232, a3 = 0.012604
2061              
2062             This routine gives the same result as the routine B in Octave 3.6.2.
2063             See also L.
2064              
2065              
2066             =cut
2067            
2068              
2069             =head2 parzen($N)
2070              
2071             The Parzen window. (Ref 1). Other names for this window are: Jackson, Valle-Poussin. This function disagrees with the Octave subroutine B, but agrees with Ref. 1.
2072             See also L.
2073              
2074              
2075             =cut
2076            
2077              
2078             =head2 parzen_octave($N)
2079              
2080             The Parzen window. No periodic version of this window is defined.
2081             This routine gives the same result as the routine B in Octave 3.6.2.
2082             See also L.
2083              
2084              
2085             =cut
2086            
2087              
2088             =head2 poisson($N,$alpha)
2089              
2090             The Poisson window. (Ref 1). This window is defined by
2091              
2092             exp (-$alpha * abs arr),
2093              
2094             where the points in arr range from -1 through 1.
2095              
2096              
2097             =cut
2098            
2099              
2100             =head2 rectangular($N)
2101              
2102             The Rectangular window. (Ref 1). Other names for this window are: dirichlet, boxcar.
2103              
2104             =cut
2105            
2106              
2107             =head2 triangular($N)
2108              
2109             The Triangular window. This window is defined by
2110              
2111             1 - abs arr,
2112              
2113             where the points in arr range from -$N/($N-1) through $N/($N-1).
2114             See also L.
2115              
2116              
2117             =cut
2118            
2119              
2120             =head2 tukey($N,$alpha)
2121              
2122             The Tukey window. (Ref 1). Another name for this window is the tapered cosine window.
2123              
2124             =cut
2125            
2126              
2127             =head2 welch($N)
2128              
2129             The Welch window. (Ref 1). Other names for this window are: Riez, Bochner, Parzen, parabolic. This window is defined by
2130              
2131             1 - arr**2,
2132              
2133             where the points in arr range from -1 through 1.
2134              
2135              
2136             =cut
2137            
2138              
2139             # Maxima code to convert between powers of cos and multiple angles in cos
2140             #grind(trigsimp(trigexpand(a0 - a1*cos(x) +a2*cos(2*x) -a3*cos(3*x) + a4*cos(4*x) -a5*cos(5*x) +a6*cos(6*x))));
2141             #
2142             #32*a6*cos(x)^6-16*a5*cos(x)^5+(8*a4-48*a6)*cos(x)^4+(20*a5-4*a3)*cos(x)^3
2143             # +(18*a6-8*a4+2*a2)*cos(x)^2+(-5*a5+3*a3-a1)*cos(x)-a6+a4-a2+a0$
2144              
2145             #(%i37) grind(trigsimp(trigreduce(c0 + c1*cos(x)
2146             # + c2*cos(x)^2 + c3*cos(x)^3 + c4*cos(x)^4 + c5*cos(x)^5 + c6*cos(x)^6)));
2147             #
2148             #(c6*cos(6*x)+2*c5*cos(5*x)+(6*c6+4*c4)*cos(4*x)+(10*c5+8*c3)*cos(3*x)
2149             # +(15*c6+16*c4+16*c2)*cos(2*x)+(20*c5+24*c3+32*c1)*cos(x)+10*c6+12*c4+16*c2+32*c0) /32$
2150              
2151             =head1 AUXILIARY SUBROUTINES
2152              
2153             These subroutines are used internally, but are also available for export.
2154              
2155             =head2 cos_mult_to_pow
2156              
2157             Convert Blackman-Harris coefficients. The BH windows are usually defined via coefficients
2158             for cosines of integer multiples of an argument. The same windows may be written instead
2159             as terms of powers of cosines of the same argument. These may be computed faster as they
2160             replace evaluation of cosines with multiplications.
2161             This subroutine is used internally to implement the Blackman-Harris
2162             family of windows more efficiently.
2163              
2164             This subroutine takes between 1 and 7 numeric arguments a0, a1, ...
2165              
2166             It converts the coefficients of this
2167              
2168             a0 - a1 cos(arg) + a2 cos( 2 * arg) - a3 cos( 3 * arg) + ...
2169              
2170             To the cofficients of this
2171              
2172             c0 + c1 cos(arg) + c2 cos(arg)**2 + c3 cos(arg)**3 + ...
2173              
2174             =head2 cos_pow_to_mult
2175              
2176             This function is the inverse of L.
2177              
2178             This subroutine takes between 1 and 7 numeric arguments c0, c1, ...
2179              
2180             It converts the coefficients of this
2181              
2182             c0 + c1 cos(arg) + c2 cos(arg)**2 + c3 cos(arg)**3 + ...
2183              
2184             To the cofficients of this
2185              
2186             a0 - a1 cos(arg) + a2 cos( 2 * arg) - a3 cos( 3 * arg) + ...
2187              
2188             =cut
2189              
2190             sub cos_pow_to_mult {
2191 0     0 1 0 my( @cin ) = @_;
2192 0 0       0 barf "cos_pow_to_mult: number of args not less than 8." if @cin > 7;
2193 0         0 my $ex = 7 - @cin;
2194 0         0 my @c = (@cin, (0) x $ex);
2195 0         0 my (@as) = (
2196             10*$c[6]+12*$c[4]+16*$c[2]+32*$c[0], 20*$c[5]+24*$c[3]+32*$c[1],
2197             15*$c[6]+16*$c[4]+16*$c[2], 10*$c[5]+8*$c[3], 6*$c[6]+4*$c[4], 2*$c[5], $c[6]);
2198 0         0 foreach (1..$ex) {pop (@as)}
  0         0  
2199 0         0 my $sign = -1;
2200 0         0 foreach (@as) { $_ /= (-$sign*32); $sign *= -1 }
  0         0  
  0         0  
2201 0         0 @as;
2202             }
2203              
2204             =head2 chebpoly
2205              
2206             =for usage
2207              
2208             chebpoly($n,$x)
2209              
2210             =for ref
2211              
2212             Returns the value of the C<$n>-th order Chebyshev polynomial at point C<$x>.
2213             $n and $x may be scalar numbers, pdl's, or array refs. However,
2214             at least one of $n and $x must be a scalar number.
2215              
2216             All mixtures of pdls and scalars could be handled much more
2217             easily as a PP routine. But, at this point PDL::DSP::Windows
2218             is pure perl/pdl, requiring no C/Fortran compiler.
2219              
2220             =cut
2221              
2222             sub chebpoly {
2223 5 50   5 1 2410 barf 'chebpoly: Two arguments expected. Got ' .scalar(@_) ."\n" unless @_==2;
2224 5         11 my ($n,$x) = @_;
2225 5 100       13 if (ref($x)) {
2226 4         14 $x = topdl($x);
2227 4 50       69 barf "chebpoly: neither $n nor $x is a scalar number" if ref($n);
2228 4         13 my $tn = zeroes($x);
2229 4         207 my ($ind1,$ind2);
2230 4         12 ($ind1,$ind2) = which_both(abs($x) <= 1);
2231 4         384 $tn->index($ind1) .= cos($n*(acos($x->index($ind1))));
2232 4         285 $tn->index($ind2) .= cosh($n*(acosh($x->index($ind2))));
2233 4         197 return $tn;
2234             }
2235             else {
2236 1 50       5 $n = topdl($n) if ref($n);
2237 1 50       5 return cos($n*(acos($x))) if abs($x) <= 1;
2238 1         13 return cosh($n*(acosh($x)));
2239             }
2240             }
2241              
2242              
2243             sub cos_mult_to_pow {
2244 0     0 1   my( @ain ) = @_;
2245 0 0         barf("cos_mult_to_pow: number of args not less than 8.") if @ain > 7;
2246 0           my $ex = 7 - @ain;
2247 0           my @a = (@ain, (0) x $ex);
2248 0           my (@cs) = (
2249             -$a[6]+$a[4]-$a[2]+$a[0], -5*$a[5]+3*$a[3]-$a[1], 18*$a[6]-8*$a[4]+2*$a[2], 20*$a[5]-4*$a[3],
2250             8*$a[4]-48*$a[6], -16*$a[5], 32*$a[6]);
2251 0           foreach (1..$ex) {pop (@cs)}
  0            
2252 0           @cs;
2253             }
2254              
2255             =head1 REFERENCES
2256              
2257             =over
2258              
2259             =item 1
2260              
2261             Harris, F.J. C,
2262             I, 1978, vol 66, pp 51-83.
2263              
2264             =item 2
2265              
2266             Nuttall, A.H. C, I,
2267             1981, vol. ASSP-29, pp. 84-91.
2268              
2269             =back
2270              
2271             =head1 AUTHOR
2272              
2273             John Lapeyre, C<< >>
2274              
2275             =head1 ACKNOWLEDGMENTS
2276              
2277             For study and comparison, the author used documents or output from:
2278             Thomas Cokelaer's spectral analysis software; Julius O Smith III's
2279             Spectral Audio Signal Processing web pages; AndrĂ© Carezia's
2280             chebwin.m Octave code; Other code in the Octave signal package.
2281              
2282             =head1 LICENSE AND COPYRIGHT
2283              
2284             Copyright 2012 John Lapeyre.
2285              
2286             This program is free software; you can redistribute it and/or modify it
2287             under the terms of either: the GNU General Public License as published
2288             by the Free Software Foundation; or the Artistic License.
2289              
2290             See http://dev.perl.org/licenses/ for more information.
2291              
2292             This software is neither licensed nor distributed by The MathWorks, Inc.,
2293             maker and liscensor of MATLAB.
2294              
2295             =cut
2296              
2297             1; # End of PDL::DSP::Windows.pm\n";
2298