Practical Least Squares FIR filter design: How to find a realistic desired response I'm using Matt. L.'s lslevin.m function to design long ($N\approx 400$) filters to an arbitrary phase and magnitude specification. The specification is the result of a numeric optimization and does cover frequencies between 200Hz and 19kHz. I found that the designed filters show poor approximation of my desired frequency response $H_1[k]$:



Apparently my desired response it hard to meet. Making the filters longer only amplifies the ripple. So I added a linear phase term $H_2[k] = H_1[k]e^{-j\frac{N-1}{2}\omega[k]}$ in the hope to make the desired response a little more realistic to meet. I am not sure why this is requried or why it makes things better, but it appears to be doing something. (Adding a (predictable) phase term to my filters is okay in my scenario as long as two filters designed with the same length have the same phase term added) Here's the result:



However, outside the range of specified frequencies I see huge spikes in the magnitude response (which is to be expected, because I don't control those regions at all):



So I extended my specification by padding $H_2[k]$ with zeros on both sides. I give a weight of 2 to all points inside my specification and a weight of 1 to all the padded zeros - my hope is that I will get a band-pass behaviour with my desired phase and magnitude response in the passband and some sort of suppression outside the passband. I don#t really care about those regions but the huge magnitude from the previous example makes the impulse response have very large coefficients with varying sign which is impossible to implement. With the padding added, things start to look better.



As you can (maybe) see, there is some ripple above 19kHz but in this case it stays low. If I increase the filter length $N$ to 600 or more - hoping to increase the accuracy of LS approximation in the passband - those ripples become extremely large, with the maxima reaching up to +50dB. Which makes this filter impossible to implement.

So as you can tell, I'm not entirely sure what I'm doing. How can I find a gradient of the linear phase term that yields the best approximation? How can I prevent huge magnitudes outside my important region of 200Hz-19kHz?
In general you seem to be doing all the right things.

The linear phase term is required since an arbitrary magnitude and phase response is typically non-causal. That means it's impossible to fit with a causal filter and your fitting error is very high. 

One way to estimate a "good" delay value is to sample your target on an FFT grid that's relatively dense and do an inverse FFT. Ideally this looks like an actual impulse response with a clear peak someplace and trailing off at the edge. Circulate this until you've captured the bulk of the energy. 

The other options is to try brute force: just try a lot of delays and pick the one that gives you the smallest error. 

Controlling the out of band behavior you can play around with the weights of the 
"don't care" bands. It's probably best to reduce the weights until the out of band energy is just about the same as the in-band energy. If you want real bandpass behavior, just apply a regular bandpass on top of your least square design. Trying to get a decent high pass at 250 Hz will cost you a lot of coefficients otherwise.

Your target has a lot of detail at low frequencies which determines the minimum number of FIR coefficients. You could try to prefit the major low frequency wiggles with an IIR filter and then just use an FIR for the remaining difference to the target.