We should add a function for computing the cross spectral density and magnitude squared coherence. Or maybe we could just bring @simonster's [Synchrony.jl]https://github.com/simonster/Synchrony.jl) implementations to DSP.jl?

High order IIR filters can be numerically unstable, therefore they are usually converted to cascades of second-order sections. The answer for this question in the DSP Stack Exchange gives a good idea on how to implement this.

freqz(digitalfilter(Lowpass(0.5), Butterworth(12)), 0)

returns 1.0 + 0.0im , but I believe it should return 0.0 + 1.0im. (Checked with Scipy)

the julia package listing currently only provides the package name and description, which for many packages, including DSP ("a DSP module for Julia"), are redundant. we should at least spell out "digital signal processing" in the description, to be more, well, descriptive! and delete "module" and "julia" since that is obvious. i would suggest squeezing "periodogram" and "window function" and "filter design" in there too, but hopefully one day the package listing page will provide a search over the readme / docs.

Some of the window functions return 2d Arrays (like rect) while others return 1d Arrays (like hanning). The issue is with using ones(n,1) and zeros(n,1) but not ones(n) and zeros(n)

@jfsantos, can you give me owner access to JuliaDSP or flip the Travis switch on this repo, and I'll add a .travis.yml file?

Now that we have @JayKickliter's excellent multirate work, it would be great to have a high-performance resampling function. I think this works to design an appropriate filter:

function resample_filter(ratio::Rational{Int}, transitionwidth, attenuation)
    # Determine resampling filter order
    L, beta = kaiserord(transitionwidth/num(ratio), attenuation)

    # Set cutoff so that we will attain `attenuation` dB rolloff
    cutoff = 1/max(num(ratio), den(ratio))-2*transitionwidth/num(ratio)

    # Design filter
    h = digitalfilter(Lowpass(cutoff), FIRWindow(kaiser(L, beta)))
    scale!(h, num(ratio))
end

The tricky part appears to be applying the filter such that the first output sample after ramp up corresponds to the first input sample (and the center of the filter).

i don't see wavelets, or principal/independent component analysis anywhere in the julia codebase or packages. did i miss them, or should we think about adding them to DSP at some point?

At present we design filters as:

digitalfilter(Lowpass(freq), Butterworth(order))
digitalfilter(Bandpass(freq1, freq2), Butterworth(order))

or equivalent, i.e., we specify filters in terms of the cutoff frequency, order, and optional filter parameters (for Chebyshev and Elliptic filters). It would be nice if we could also specify filters in terms of the cutoff frequency, stopband frequency, passband ripple, and stopband attenuation. We could have something like:

digitalfilter(Lowpass(0.2; Rp=1, Ws=0.5, Rs=90), Butterworth())
digitalfilter(Highpass(0.5; Rp=1, Ws=0.2, Rs=90), Butterworth())

Then for bandpass filters we'd pass Ws and (optionally) Rs as tuples:

digitalfilter(Bandpass(0.2, 0.3; Rp=1,
                       Ws=(0.1, 0.5), Rs=(60, 90)),
              Butterworth())

And for bandstop filters, we'd specify the passband frequencies Wp instead of Ws:

digitalfilter(Bandstop(0.2, 0.3; Rp=(.5, 1),
                       Wp=(0.1, 0.4), Rs=90),
              Butterworth())

This seems a bit complex, but I haven't thought of anything better yet.

the normalization "so that the area under the periodogram is equal to the uncentered variance of the original signal" is rather unusual.
Possibility to have also the classical normalisation (f_e/N)*abs(FFT)^2 would be great !

Here's a rough draft. I have next to zero knowledge of linear algebra, or experience using this design method. Filing this as an issue for now, in case someone with more experience wants to take a crack at it.

function firls(n::Integer, bands::AbstractVector)
    isodd(n) || throw(ArgumentError("Number of taps must be odd"))

    m = int((n-1)/2)
    Q = zeros(m+1,m+1)
    b = zeros(m+1,m+1)

    for (ƒl, ƒh, amplitude) in bands
        for j in 0:m 

            b[j+1] += amplitude*2π*(ƒh*sinc(ƒh*j) - ƒl*sinc(ƒl*j))

            for i in 0:m
                Q[i+1,j+1] += π*( ƒh*sinc(ƒh*(i-j)) - ƒl*sinc(ƒl*(i-j)) ) # Q₁, Toeplitz
                Q[i+1,j+1] += π*( ƒh*sinc(ƒh*(i+j)) - ƒl*sinc(ƒl*(i+j)) ) # Q₂, Hankel
            end
        end
    end

    a = Q \ b
    h = [a[m+1:-1:2]./2, a[1], a[2:m+1]./2]
end

# Lowpass
N    = 41
b    = [(0.0, 0.4, 1), (0.6, 1, 0)]
h_lp = firls(N, b)

# Highpass
N    = 41
b    = [(0.0, 0.4, 0), (0.6, 1, 1)]
h_hp = firls(N, b)

# Bandpass
N    = 41
b    = [(0, .25, -50dBa), (.4, .6, 0dBa), (.75, 1, -50dBa)]
h_bp = firls(N, b)

# Bandstop
N    = 41
b    = [(0, .2, 0dBa), (.4, .6, -10dBa), (.8, 1, 0dBa)]
h_bs = firls(N, b)

Theory
SciPy discussion

Is there any interest in an equivalent matlab filtfilt function?
Or is this a bit high level and people can implement their own?

I thinks these are not included yet? I could contribute them if you wish.

We need better documentation for the functions here than a list in the README. We should probably pick some other package whose docs we like and use the same approach they do. It looks like most packages have one page per category of functions (e.g. Distributions) although some have one page per function (e.g. Winston). The Julia docs have one page for all functions in the standard library, although that approach seems less than optimal for functions that do more complicated things.

it would be nice to add checks on the values of parameters passed into dpss, much like in gaussian and tukey. i know that ntapers must be <= 2*nw-1, which is what we have as a default. matlab tells me that nw must be less than n/2, yet the julian dpss as it currently stands does not produce an error given dpss(9,5,9). if someone can confirm what the limit should be, i'd be happy to submit a PR.

Currently,

there is no functionality implemented that computes and returns an array for the frequencies of the DFT. I think there are two possibilities which I discuss here to start the discussion on the interface.

1. make periodogram, welch_pgram and bartlett_pgram return two arrays by default, one for the frequencies and one for the power.

1b) add an option that triggers these function to not return the frequency array.

2. simply add a new function that computes the required frequency array, which is called directly by the user.

Are there any functions yet to calculate a filter response at set frequencies similar to matlabs freqz and freqs?

It seems that DSP.jl and Images.jl define conflicting versions of (-) and (+):

julia> import Images

julia> import DSP
Warning: New definition 
    -(AbstractArray{T,1},Frequencies) at /home/kevin/.julia/v0.3/DSP/src/util.jl:71
is ambiguous with: 
    -(AbstractImageDirect{T,N},AbstractArray{T,N}) at /home/kevin/.julia/v0.3/Images/src/algorithms.jl:25.
To fix, define 
    -(AbstractImageDirect{T,1},Frequencies)
before the new definition.
Warning: New definition 
    +(AbstractArray{T,1},Frequencies) at /home/kevin/.julia/v0.3/DSP/src/util.jl:71
is ambiguous with: 
    +(AbstractImageDirect{T,N},AbstractArray{T,N}) at /home/kevin/.julia/v0.3/Images/src/algorithms.jl:14.
To fix, define 
    +(AbstractImageDirect{T,1},Frequencies)
before the new definition.

(Interestingly, if DSP is loaded first, only (-) conflicts...)

This was reported by @jcbsnd in #104 (comment) and #104 (comment); I'm splitting it off from that issue because it's not strictly related to resampling.

I confirmed that designing and applying the corresponding SOS filter in MATLAB gives exactly the same results as in DSP.jl, so there is no bug in our implementation, but it does seem like some parameters should be tweaked.

I tried extending the length of the extrapolation at the beginning and end of the signal so that it is 6*nbiquads instead of 6, and that substantially improved the result:

I don't really see any reason not to do this, but I guess I should also try with cases where the first sample is not 0.

Right now we test the output types of some of the windows, but we should also check that the results are correct.

using DSP (master branch)
Julia Version 0.3.0-prerelease+2954

filt! not defined
while loading /home/berceanu/.julia/v0.3/DSP/src/filter_design.jl, in expression starting on line 259
while loading /home/berceanu/.julia/v0.3/DSP/src/DSP.jl, in expression starting on line 7

http://iainnz.github.io/packages.julialang.org/?pkg=DSP&ver=0.2

We need to at least add some reference filters/filtering results so that we can compare the values computed by our functions with them. Another option is to use PyCall and compare our results with Scipy (but then we'll depend on the correctness of Scipy's code).

I'm pretty sure I remember using welch_pgram in the past with complex vectors. Is this the change in FFTW.jl Plan API I read about in a different issue?

julia> welch_pgram(rand(Complex64,100))
ERROR: `Plan{T<:Union(Float64,Complex{Float64},Complex{Float32},Float32)}` has no method matching Plan{T<:Union(Float64,Complex{Float64},Complex{Float32},Float32)}(::Array{Complex{Float32},1}, ::Array{Complex{Float32},1}, ::Int64, ::Uint32, ::Float64)
 in welch_pgram at /Users/jkickliter/.julia/v0.3/DSP/src/periodograms.jl:288
 in welch_pgram at /Users/jkickliter/.julia/v0.3/DSP/src/periodograms.jl:280 (repeats 2 times)

julia> spectrogram(rand(Complex64,100))
ERROR: `Plan{T<:Union(Float64,Complex{Float64},Complex{Float32},Float32)}` has no method matching Plan{T<:Union(Float64,Complex{Float64},Complex{Float32},Float32)}(::Array{Complex{Float32},1}, ::Array{Complex{Float32},1}, ::Int64, ::Uint32, ::Float64)
 in spectrogram at /Users/jkickliter/.julia/v0.3/DSP/src/periodograms.jl:325
 in spectrogram at /Users/jkickliter/.julia/v0.3/DSP/src/periodograms.jl:316 (repeats 2 times)

should consider adding this as an option somehow

Most implementations of spectrogram return both the STFT and the frame-wise PSD, but ours returns only the PSDs. Should we adapt our Spectrogram type to add the STFT to it as well or create a separate STFT function (which would then be used by spectrogram)?

INFO: Loading help data...
Warning: New definition
-(AbstractArray{T,1},Frequencies) at /home/berceanu/.julia/v0.3/DSP/src/util.jl:71
is ambiguous with:
-(AbstractImageDirect{T,N},AbstractArray{T,N}) at /home/berceanu/.julia/v0.3/Images/src/algorithms.jl:25.
To fix, define
-(AbstractImageDirect{T,1},Frequencies)
before the new definition.
Warning: New definition
+(AbstractArray{T,1},Frequencies) at /home/berceanu/.julia/v0.3/DSP/src/util.jl:71
is ambiguous with:
+(AbstractImageDirect{T,N},AbstractArray{T,N}) at /home/berceanu/.julia/v0.3/Images/src/algorithms.jl:14.
To fix, define
+(AbstractImageDirect{T,1},Frequencies)

Our implementation of the Hann window uses [0.5*(1 - cos(2*pi*k/(n-1))) for k=0:(n-1)], so it always returns symmetric windows. MATLAB and Scipy offer the option of using an asymmetric window, which is used in spectral analysis. Should we add another function for asymmetric windows or do something similar to what is done in Scipy? (i.e., check for a flag and whether the length of the window is odd/even, and generate the appropriate window)

This issue is being filed by a script, but if you reply, I will see it.

PackageEvaluator.jl is a script that runs nightly. It attempts to load all Julia packages and run their test (if available) on both the stable version of Julia (0.2) and the nightly build of the unstable version (0.3).

The results of this script are used to generate a package listing enhanced with testing results.

The status of this package, DSP, on...

• Julia 0.2 is 'Tests fail, but package loads.'
• Julia 0.3 is 'Tests fail, but package loads.'

'No tests, but package loads.' can be due to their being no tests (you should write some if you can!) but can also be due to PackageEvaluator not being able to find your tests. Consider adding a test/runtests.jl file.

'Package doesn't load.' is the worst-case scenario. Sometimes this arises because your package doesn't have BinDeps support, or needs something that can't be installed with BinDeps. If this is the case for your package, please file an issue and an exception can be made so your package will not be tested.

This automatically filed issue is a one-off message. Starting soon, issues will only be filed when the testing status of your package changes in a negative direction (gets worse). If you'd like to opt-out of these status-change messages, reply to this message.

can someone please explain why Base.time gets overwritten by using DSP unless Base.time is accessed beforehand? so for example, in a brand new session of julia:

julia> methods(time)
#2 methods for generic function "time":
time(tm::TmStruct) at libc.jl:70
time() at util.jl:4

julia> using DSP

julia> methods(time)
#2 methods for generic function "time":
time(tm::TmStruct) at libc.jl:70
time() at util.jl:4

julia> time()
1.408481824646962e9

but if i ctrl-D out of that session and start a new one and then:

julia> using DSP

julia> methods(time)
#1 method for generic function "time":
time(p::Spectrogram{T,F<:Union(Frequencies,Range{T})}) at /home/arthurb/.julia/v0.4/DSP/src/periodograms.jl:310

julia> time()
ERROR: `time` has no method matching time()
you may have intended to import Base.time

if i then go edit periodograms.jl to add import Base.time, then start a new session:

julia> using DSP

julia> methods(time)
#3 methods for generic function "time":
time(tm::TmStruct) at libc.jl:70
time() at util.jl:4
time(p::Spectrogram{T,F<:Union(Range{T},Frequencies)}) at /home/arthurb/.julia/v0.4/DSP/src/periodograms.jl:311

julia> time()
1.408481871189969e9

seems this is the behavior we want, and it's an easy fix, and i'll gladly submit a PR, but i don't understand why Base.time gets overwritten in the first place. the docs say that using just makes available the exported items for searching against undefined names. but time is already defined in Base, so why is it being overwritten? is this a bug in julia or a misunderstanding on my part?

I'd like to publish an initial release of this package in METADATA.jl so that people know it's out there. @jfsantos, any objections?

I want to flesh out the docs a little so that it's clearer how to actually design filters, but before I do that, I think it might make sense to reorganize the filter-related functions. Right now we have:

• FFTFilt
  • fftfilt and firfilt (should probably have mutating versions of these as well)
• FilterDesign
  • ZPKFilter, TFFilter, SOSFilter, and BiquadFilter types, conversions between types, and filt implementations
  • Butterworth filter prototype
  • Prototype transformation and bilinear transform
• FilterResponse
  • freqz and freqs
• ZeroPhaseFiltering
  • filtfilt

All of these modules but FFTFilt depend on the FilterDesign module because it defines the filter types, which suggests that this should just be one module. I propose a single DSP.Filter module that would reside in a separate subdirectory and contain the following files:

• types.jl - ZPKFilter, TFFilter, SOSFilter, and BiquadFilter type definitions. When we have FIR filter design we will probably want an FIRFilter type as well. Not sure if the conversions should live here or in a separate file.
• filt.jl - filt, fftfilt, firfilt, and filtfilt implementations (eventually firfilt will probably just be filt(::FIRFilter, x))
• design.jl - prototypes, transformations, and bilinear transform
• response.jl - freqz and freqs

Does this sound good? Any objections?

PackageEvaluator.jl is a script that runs nightly. It attempts to load all Julia packages and run their tests (if available) on both the stable version of Julia (0.3) and the nightly build of the unstable version (0.4). The results of this script are used to generate a package listing enhanced with testing results.

On Julia 0.4

• On 2014-09-29 the testing status was Tests pass.
• On 2014-09-30 the testing status changed to Tests fail, but package loads.

Tests pass. means that PackageEvaluator found the tests for your package, executed them, and they all passed.

Tests fail, but package loads. means that PackageEvaluator found the tests for your package, executed them, and they didn't pass. However, trying to load your package with using worked.

This error on Julia 0.4 is possibly due to recently merged pull request JuliaLang/julia#8493.
This issue was filed because your testing status became worse. No additional issues will be filed if your package remains in this state, and no issue will be filed if it improves. If you'd like to opt-out of these status-change messages, reply to this message saying you'd like to and @IainNZ will add an exception. If you'd like to discuss PackageEvaluator.jl please file an issue at the repository. For example, your package may be untestable on the test machine due to a dependency - an exception can be added.

Test log:

>>> 'Pkg.add("DSP")' log
INFO: Installing DSP v0.0.4
INFO: Installing Polynomials v0.0.3
INFO: Installing Reexport v0.0.1
INFO: Package database updated

>>> 'using DSP' log
Julia Version 0.4.0-dev+856
Commit 46c7bbf (2014-09-30 04:44 UTC)
Platform Info:
  System: Linux (x86_64-unknown-linux-gnu)
  CPU: Intel(R) Xeon(R) CPU E5-2650 0 @ 2.00GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas
  LIBM: libopenlibm
  LLVM: libLLVM-3.3

>>> test log

ERROR: `ZPKFilter{Z<:Number,P<:Number,K<:Number}` has no method matching ZPKFilter{Z<:Number,P<:Number,K<:Number}(::Array{Any,1}, ::Array{Complex{Float64},1}, ::Int64)
 in include at ./boot.jl:245
 in include_from_node1 at ./loading.jl:128
 in include at ./boot.jl:245
 in include_from_node1 at loading.jl:128
 in process_options at ./client.jl:285
 in _start at ./client.jl:354
 in _start_3B_3624 at /home/idunning/julia04/usr/bin/../lib/julia/sys.so
while loading /home/idunning/pkgtest/.julia/v0.4/DSP/test/filter_design.jl, in expression starting on line 10
while loading /home/idunning/pkgtest/.julia/v0.4/DSP/test/runtests.jl, in expression starting on line 4
INFO: Testing DSP
=================================[ ERROR: DSP ]=================================

failed process: Process(`/home/idunning/julia04/usr/bin/julia /home/idunning/pkgtest/.julia/v0.4/DSP/test/runtests.jl`, ProcessExited(1)) [1]

================================================================================
INFO: No packages to install, update or remove
ERROR: DSP had test errors
 in error at error.jl:21
 in test at pkg/entry.jl:719
 in anonymous at pkg/dir.jl:28
 in cd at ./file.jl:20
 in cd at pkg/dir.jl:28
 in test at pkg.jl:68
 in process_options at ./client.jl:213
 in _start at ./client.jl:354
 in _start_3B_3624 at /home/idunning/julia04/usr/bin/../lib/julia/sys.so


>>> end of log

Many functions in this library (for example, all FFT-related functions) are part of the 0.2.0 standard library. These functions have to be removed and the remaining code has to run correctly using the functions provided by the standard lib.

fftshift(fftfreq(10)) fails with

type: circshift: in typeassert, expected Frequencies, got Array{Float64,1}

E.g., if I want a 2D tukey window, I would do

win = tukey((w,h), 0.5)

which would be equivalent to

win = tukey(w, 0.5)*tukey(h, 0.5)'

These are in some sense trivial, and it's questionable whether they belong here or Images.jl, or if as a user I should just define them on my own. In truth, I'm not even sure how standard these are, but I was finding myself in need of the filter above.

Anyway, I figured I'd start here and get feedback.

Cc: @timholy

Useful notes here. I implemented Chebyshev type I and II filter design locally, but if we implement elliptic filter design, it seems like we could just treat these as special cases.

Add code for adaptive (LMS/RLS) filters. In case this requires adding too much to this package, it could be done in a separate package instead.

PackageEvaluator.jl is a script that runs nightly. It attempts to load all Julia packages and run their tests (if available) on both the stable version of Julia (0.3) and the nightly build of the unstable version (0.4). The results of this script are used to generate a package listing enhanced with testing results.

On Julia 0.4

• On 2015-01-22 the testing status was Tests pass.
• On 2015-01-23 the testing status changed to Tests fail, but package loads.

Tests pass. means that PackageEvaluator found the tests for your package, executed them, and they all passed.

Tests fail, but package loads. means that PackageEvaluator found the tests for your package, executed them, and they didn't pass. However, trying to load your package with using worked.

This issue was filed because your testing status became worse. No additional issues will be filed if your package remains in this state, and no issue will be filed if it improves. If you'd like to opt-out of these status-change messages, reply to this message saying you'd like to and @IainNZ will add an exception. If you'd like to discuss PackageEvaluator.jl please file an issue at the repository. For example, your package may be untestable on the test machine due to a dependency - an exception can be added.

Test log:

>>> 'Pkg.add("DSP")' log
INFO: Installing DSP v0.0.6
INFO: Installing Polynomials v0.0.3
INFO: Installing Reexport v0.0.2
INFO: Package database updated

>>> 'using DSP' log
Julia Version 0.4.0-dev+2876
Commit f164ac1 (2015-01-22 22:58 UTC)
Platform Info:
  System: Linux (x86_64-unknown-linux-gnu)
  CPU: Intel(R) Xeon(R) CPU E5-2650 0 @ 2.00GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas
  LIBM: libopenlibm
  LLVM: libLLVM-3.3

>>> test log

WARNING: iceil(x) is deprecated, use ceil(Integer,x) instead.
 in depwarn at ./deprecated.jl:40
 in iceil at deprecated.jl:29
 in dpss at /home/idunning/pkgtest/.julia/v0.4/DSP/src/windows.jl:109
 in include at ./boot.jl:249
 in include_from_node1 at ./loading.jl:128
 in include at ./boot.jl:249
 in include_from_node1 at loading.jl:128
 in process_options at ./client.jl:319
 in _start at ./client.jl:403
WARNING: iceil(x) is deprecated, use ceil(Integer,x) instead.
 in depwarn at ./deprecated.jl:40
 in iceil at deprecated.jl:29
 in dpss at /home/idunning/pkgtest/.julia/v0.4/DSP/src/windows.jl:109
 in anonymous at test.jl:85
 in do_test at test.jl:47
 in include at ./boot.jl:249
 in include_from_node1 at ./loading.jl:128
 in include at ./boot.jl:249
 in include_from_node1 at loading.jl:128
 in process_options at ./client.jl:319
 in _start at ./client.jl:403
ERROR: LoadError: LoadError: DimensionMismatch("first array argument must be the same length as the second array argument")
 in generic_scale! at linalg/generic.jl:24
 in iir_filtfilt at /home/idunning/pkgtest/.julia/v0.4/DSP/src/Filters/filt.jl:122
 in filtfilt at /home/idunning/pkgtest/.julia/v0.4/DSP/src/Filters/filt.jl:168
 in include at ./boot.jl:249
 in include_from_node1 at ./loading.jl:128
 in include at ./boot.jl:249
 in include_from_node1 at loading.jl:128
 in process_options at ./client.jl:319
 in _start at ./client.jl:403
while loading /home/idunning/pkgtest/.julia/v0.4/DSP/test/filt.jl, in expression starting on line 167
while loading /home/idunning/pkgtest/.julia/v0.4/DSP/test/runtests.jl, in expression starting on line 6

INFO: Testing DSP
=================================[ ERROR: DSP ]=================================

failed process: Process(`/home/idunning/julia04/usr/bin/julia --check-bounds=yes --code-coverage=none --color=no /home/idunning/pkgtest/.julia/v0.4/DSP/test/runtests.jl`, ProcessExited(1)) [1]

================================================================================
INFO: No packages to install, update or remove
ERROR: DSP had test errors
 in error at error.jl:19
 in test at pkg/entry.jl:717
 in anonymous at pkg/dir.jl:28
 in cd at ./file.jl:20
 in cd at pkg/dir.jl:28
 in test at pkg.jl:68
 in process_options at ./client.jl:241
 in _start at ./client.jl:403

>>> end of log

Currently our ArraySplit type does not zero-pad the signal anywhere, but one may wish to zero pad the beginning or the end of an array (for example, if computing MDCTs or an inverse STFT). Should we add this functionality to ArraySplit?

A straightforward way of adding it would be to have two attributes in ArraySplit that store the number of zeros at the beginning and the end, and in ArraySplit.getindex check if the required range is in one of these intervals (in which case we would just fill the corresponding elements with zeros). This may be bad performance-wise, but we can try to avoid having to do a comparison prior to copying each element to the "view" buffer.

should consider adding this as an option somehow.

Add FIR filter design functions (window method, Remez exchange algorithm). These functions could be ported from Scipy.

The documentation here appears to be using s=z^-1, which I haven't seen before. In signal processing and controls, I have always seen s to mean the Laplace-transform domain variable. Perhaps something like D=z^-1 can be used, where D is called the "delay operator".

It might also be confusing that z is used as the roots of the numerator polynomial of a transfer function and also the z-transform domain variable as in H(z)

I made a few tweeks to spectrogram which results in quite significant 4 times faster method using half the memory. Also, by changing the window to w=(n)->ones(n) the output is a 1d Array if the input is a 1d Array.

EDIT: I got 40 times speedup and used 5th of the memory:

spectrogramf
elapsed time: 0.017406119 seconds (6052416 bytes allocated)
Array{Float64,1}
4096
spectrogramf2
elapsed time: 0.002346886 seconds (2663856 bytes allocated)
Array{Float64,1}
4096
spectrogram
elapsed time: 0.084957317 seconds (11150800 bytes allocated)
Array{Float64,2}
4096
welch_pgram
elapsed time: 0.040620592 seconds (7973488 bytes allocated, 55.53% gc time)
Array{Float64,1}
4096

using DSP, Base.Test
using Base.FFTW: fftwNumber
typealias dspNumber fftwNumber

function spectrogramf{T<:dspNumber}(s::Array{T,1}; n=int(length(s)/8), m=int(n/2), r=1, w=(n)->ones(n))
  w = w(n)
  ind = arrayspliti(s, n, m)
  sf = zeros(eltype(s),n)
  for i=ind
    sf += periodogramf(s[i:i+n-1].*w)
  end

  t=( (0:length(ind)-1)*(n-m) + n/2) / r
  f=(0:n-1)/n*r
  sf/(r*length(ind)), t, f
end
function spectrogramf2{T<:dspNumber}(s::Array{T,1}; n=int(length(s)/8), m=int(n/2), r=1, w=(n)->ones(n))
  w = w(n)
  ind = arrayspliti(s, n, m)
  sf = zeros(eltype(s),n)
  per = zeros(typeof(1.0im),n)
  for i=ind
    per[:] = s[i:i+n-1].*w
    periodogramtl!(per)
    sf += real(per)
  end

  t=( (0:length(ind)-1)*(n-m) + n/2) / r
  f=(0:n-1)/n*r
  sf/(r*length(ind)*length(per)), t, f
end
function periodogramf{T<:dspNumber}(s::Array{T,1})
    s_fft = fft(s)
    abs2(s_fft)/length(s_fft)
end
function periodogramtl!{T<:dspNumber}(s::Array{T,1})
    fft!(s)
    s[:] = abs2(s)
    return nothing
end

function arraysplit_su(sn::Integer, n::Integer, m::Integer)
    # n = m is a problem - the algorithm will not terminate.
    if !(0 <= m < n)
        error("m must be between zero and n.")
    end

    # the length of the non-overlapping array stride
    l = n - m

    # total number of strides is the total length of the signal divided
    # by the unique number of elements per stride. extra elements at the
    # end of of the signal are dropped.
    k = ifloor(sn/l - n/l + 1)
    l, k
end
function arrayspliti{T<:dspNumber}(s::Array{T,1}, n::Integer, m::Integer)
    l, k = arraysplit_su(length(s), n, m)
    [a*l + 1 for a=0:(k-1)]
end


# test

N=1024*32
x0=randn(N)
nn=int(N/8)
mm=int(N/16)

# initialize variables before timing
sf, tf, ff = spectrogramf(x0, n=nn, r=10, m=mm)
sf2, tf2, ff2 = spectrogramf2(x0, n=nn, r=10, m=mm)
p, t, f = spectrogram(x0, n=nn, r=10, m=mm)
sw= welch_pgram(x0, nn, mm)

println(spectrogramf)
@time begin sf, tf, ff = spectrogramf(x0, n=nn, r=10, m=mm)
end
println(typeof(sf))
println(length(sf))

println(spectrogramf2)
@time begin sf2, tf2, ff2 = spectrogramf2(x0, n=nn, r=10, m=mm)
end
println(typeof(sf))
println(length(sf))

println(spectrogram)
@time begin p, t, f = spectrogram(x0, n=nn, r=10, m=mm)
s = 1/size(p,2)*sum(p,2)
end
println(typeof(s))
println(length(s))

println(welch_pgram)
@time begin sw= welch_pgram(x0, nn, mm)
end
println(typeof(sw))
println(length(sw))


@test_approx_eq s sf2
@test_approx_eq s sf
@test_approx_eq t tf
@test_approx_eq f ff
@test_approx_eq s sw/10

Hi, it would be nice to have a function to generate a FIR filter with an arbitrary frequency response as in the following implementations:

scipy (firwin2 function):
http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.firwin2.html#scipy.signal.firwin2

octave (fir2 function in the octave-forge signal package):
http://octave.sourceforge.net/signal/function/fir2.html

R (fir2 function in the signal package):
http://cran.r-project.org/web/packages/signal/signal.pdf

With the most recent DSP and Julia master branches, DSP/test/filter_design.jl fails.

Here's the output:

$ ./julia ~/.julia/DSP/test/filter_design.jl
ERROR: no method convert(Type{Array{Complex{Float64},1}}, Array{Float64,1})
in convert at /home/rick/.julia/DSP/src/filter_desig
n.jl:109
in include_from_node1 at loading.jl:120
while loading /home/rick/.julia/DSP/test/filter_design.jl, in expression starting on line 82

$ ./julia -e 'versioninfo()'
Julia Version 0.3.0-prerelease+1249
Commit aefde34* (2014-01-26 15:48 UTC)
Platform Info:
System: Linux (i686-redhat-linux)
CPU: Genuine Intel(R) CPU T2250 @ 1.73GH
z
WORD_SIZE: 32
BLAS: libopenblas (DYNAMIC_ARCH NO_AFFINITY)
LAPACK: libopenblas
LIBM: libopenlibm

I think they confuse more than they help. The documentation for the Julia standard library and other packages seem to omit them.