Julius Smith's background is in electrical engineering (BS Rice 1975, PhD Stanford 1983). He is presently Professor of Music and (by courtesy) of Electrical Engineering at Stanford's Center for Computer Research in Music and Acoustics (CCRMA), teaching courses and pursuing research related to signal processing applied to music and audio systems. See http://ccrma.stanford.edu/~jos/ for details.

Fitting Filters to Measured Amplitude Response Data Using invfreqz in Matlab

Julius Orion Smith III October 11, 20102 comments

This blog post has been moved to the code snippet section and can now be found HERE.  Please update your bookmark.  Thanks!


Music/Audio Signal Processing

Julius Orion Smith III September 5, 20087 comments

Greetings,

This is my blog from the point of view of a music/audio DSP research engineer / educator. It is informal and largely nontechnical because nearly everything I have to say about signal processing is (or will be) somewhere in my four-book series: Mathematics of DFT with Audio Applications, Introduction to Digital Filters, Physical Audio Signal Processing and


Re: Narrow band filter at high rate

Reply posted 4 years ago (08/26/2017)
Ok, that makes sense.  To answer the question, we need to know the pass-band and stop-band edge frequencies for the desired isolation filter.  A 6us delay spec...

Re: Narrow band filter at high rate

Reply posted 4 years ago (08/26/2017)
This sounds dubious to me.  One cycle at 20 KHz is already 50 us.  Where is the 20 KHz band?  If centered on 2.5 MHz, then the "Q" of the band is Q = 2.5e6/20e3...

Re: Typical filter orders for DSPs

Reply posted 4 years ago (08/04/2017)
For custom FIR filter designs, I like to start with a gigantic, overkilled, power-of-2-length design given simply by the IFFT of the uniformly sampled desired...
I can recommend this: https://cnx.org/contents/feIAPEwX@1.1:vvZSF2za@2/S...
You should provide a small complete test example instead of the incomplete real thing.16002000 should probably be 1600 2000

Re: Magnitude of frequency components

Reply posted 4 years ago (06/28/2017)
I would not introduce a window here because the current approach is going for exactly one period of 25 Hz in the FFT input buffer.  This period should not be distorted...

Re: Estimating SNR without windowing

Reply posted 4 years ago (06/17/2017)
I would replace A by A/2 in your formula and sum over positive frequencies only (which maybe you're already doing).  The FFT bin magnitude is N*A/2 not N*A in...
I would use notch filters for this. Your largest peak looks like it might require an adaptive notch filter.  Search on these terms for details.

Re: Simplification of an LMS implementation

Reply posted 4 years ago (06/12/2017)
...

Re: Power spectral density

Reply posted 4 years ago (05/22/2017)
I would definitely use Welch's method, which is based on the periodogram.  Since your only issue is doing it sequentially through a file instead of on one giant...
Another thing to watch out for is clipping, for example if N*A*N clips, but N*A doesn't.

Re: Filtering of amplitude spectrum

Reply posted 4 years ago (04/21/2017)
For this I would use a "zero-phase" smoother (such as a moving average symmetric around its midpoint) on the COMPLEX spectrum, and then form the magnitude.Note that...
Yes, Matlab's hilbert() function expects a complete period of a periodic input signal, in this case it needs a complete period of the log-magnitude response...
I assume by "dynamic range" you mean "maximum gain".  I would use max(abs(fft(impulseResponse,veryLargeFFTSize)))Sampling-rate converters are linear, time-varying...

Re: Is there any way to prevent Deconvolution?

Reply posted 5 years ago (04/07/2017)
Lot's of 0s around the unit circle in the transfer function :-)Non-minimum-phase 0s make it a bit harder, but you can still simply Fourier transform, divide by...

Re: Real or complex FFT on IQ data

Reply posted 5 years ago (03/30/2017)
The information is the same.  Since the FFT is a linear operator, FFT(I + j*Q) = FFT(I) + j*FFT(Q).  Now break down each FFT as FFT(x) = DCT(x) - j*DST(x), and...

Re: The adoption of the frequency measure "Hertz"

Reply posted 5 years ago (03/29/2017)
I am glad to learn that the IEEE was simply following the lead of major standards bodies.  I've always thought the standard should be 'cps', and if we're going...

Re: Complex envelope and analytic signals.

Reply posted 5 years ago (03/26/2017)
For a more detailed discussion, you might tryhttps://www.dsprelated.com/freebooks/sasp/Hilbert_...

Re: Frequency resolution

Reply posted 5 years ago (03/24/2017)
Apologies for likely duplication, but the optimal window for frequency-resolution in stationary noise is the rectangular window.  See the classic papers by Rife...

Re: If H(ejw) is an ideal LPF, what is H(-ejw) ?

Reply posted 5 years ago (03/21/2017)
Negating z is a pi-rotation, so you get a high-pass filter.
I would also plot W, and limit its variation using the threshold I mentioned.
Maybe change it to W = 1./(1.0e-4 + abs(Sxy));
I would be very worried about this line:W = 1./abs(Sxy)
Generalized cross-correlation does not pertain to cross-correlating two sinusoids.  However, it's a good acid-test of the implementation which should not blow...
Check out cfirpm: https://www.mathworks.com/help/signal/ug/fir-filte...
Yes, only the noise should be whitened.  In audio, we typically use segments of silence to form a noise-floor estimate.  Otherwise, you can try to preserve the...
Many detection and estimation formulations assume that the added noise is white (see "minimum variance" and "maximum likelihood" estimation).  When the noise...

Re: harmonic distortion of a single-tone excitation

Reply posted 5 years ago (02/08/2017)
The amplitudes of the harmonic-distortion components are smaller than the fundamental amplitude for a "weak nonlinearity", but not in general.  Consider the squaring...

Re: Frequency multiplication in the digital domain

Reply posted 5 years ago (02/05/2017)
For frequency shifting, filter out the negative-frequency component and multiply by a similar "complex sinusoid" tuned to the "difference frequency", thereby summing...

Re: Frequency multiplication in the digital domain

Reply posted 5 years ago (02/05/2017)
If you know the frequency, it is easy to delay a sinusoid by a quarter cycle (and negate) to turn sin into cos.  Just make a first-order allpass filter having a...
This is easy if there is no noise, and difficult to impossible in the presence of noise.  Your "difference frequency" is 0.01 Hz, which needs 100 seconds for a...

Re: How to create a vector graphic of a block model

Reply posted 5 years ago (01/30/2017)
Yes, tikz is the gold standard.  I use xfig myself.

Re: Spurious when Fs/F is not an integer

Reply posted 5 years ago (01/23/2017)
Another common practice is to do linear interpolation based on the low-order wavetable address bits.  Below, for example, is the interpolating sinusoid oscillator...

Re: Upsampled input to an Adaptive filter?

Reply posted 5 years ago (01/23/2017)
Did the second sinc correspond to the lowpass filter or the upsampling factor?  I.e., did the zero-crossings move?  If the sinc function "stretched" horizontally,...
This is not a classical control answer, but I like to think in terms of quality factor, or Q. The Q is defined as center frequency over bandwith.  It is approximately...
In my experience, one finds the eigenvalues of the state-transition matrix (in the state-space representation) and see that they all have magnitude less than 1.

Re: Design a filter-nonlinearity model of an amplifier

Reply posted 5 years ago (01/06/2017)
This can be called a Wiener Model, for which there is support in Matlab: https://www.mathworks.com/help/ident/ug/identifying-hammerstein-wiener-models.html

Re: Audio Beamforming

Reply posted 5 years ago (11/13/2016)
Instead of switching, you can "cross-fade" over some time interval, such as 100ms.  Linear amplitude ramps are pretty standard, but I would probably use half a...
Yes, the inverse DFT can be interpreted as the inverse Fourier series (or "additive synthesis") generating the "periodic extension" of the original finite-time...

Re: A New Reverse IIR Filtering Algorithm

Reply posted 5 years ago (10/21/2016)
I've never seen that - interesting - complexity nicely between TIIR and plain FIR

Re: Proof of Z-transform of M-Fold Decimator/Compressor

Reply posted 5 years ago (10/21/2016)
The "change of variable" approach you propose can be understood from the Fourier "scaling theorem".  Instead of sampling the time domain, it will scale the time...

Re: Zero-Padding as scalloping loss attenuator

Reply posted 5 years ago (10/21/2016)
@Tim: Scalloping loss refers to the "droop" in magnitude that occurs between FFT bins when you don't combine them properly.  It depends on the window used.  For...

Re: Zero-Padding as scalloping loss attenuator

Reply posted 5 years ago (10/20/2016)
There should be no such thing as "scalloping loss" :-)Seriously, when a sinusoid's frequency falls between FFT bins, you should combine adjacent bins to get a decent estimate...

Re: Steepest-slope 2nd-order resonant filter?

Reply posted 5 years ago (10/15/2016)
How about just measuring it to see what's going on?: https://ccrma.stanford.edu/realsimple/imp_meas/I agree that a "second-order bandpass with an extremely steep...

Re: SSB Demodulation

Reply posted 5 years ago (10/11/2016)
I think you are correct.  However, by converting the imaginary part so that it can be combined with the real part, you use all of the signal power instead of...

Re: Removal of DC Spike from quadrature receiver

Reply posted 5 years ago (08/14/2016)
If you use a rectangular window encompassing the entire length of the signal, then the dc component (bin 0) of the de-meaned signal's FFT will be exactly zero,...

Re: Removal of DC Spike from quadrature receiver

Reply posted 5 years ago (08/13/2016)
Since subtracting the mean does not work, it sounds like you have more low-frequency noise than just at dc.  How about a high-pass filter having zeros at dc...
Yes, I mean simply that we need H(j 2 pi f) = 0 for all |f| greater than or equal to fs/2 (half the sampling rate = Nyquist limit).  Otherwise, we are violating...
I agree with Tim - there should be no amplitude at the Nyquist limit when simulating an analog system.  As you have discovered, it is impossible to specify both the...

Re: PLL in presence of noise

Reply posted 5 years ago (07/26/2016)
I would try phase-unwrapping (unwrap() in matlab) followed by linear regression (e.g., polyfit(time,phases,1)).  Plot the unwrapped phase to make sure it generally...

Re: "System Identification" Ideas

Reply posted 5 years ago (07/25/2016)
Sounds like a fun problem.  You should be able to apply your solution to automatic buying/selling decisions in the stock market :-)It sounds like simple characterizations...

Re: Similarity between two data points

Reply posted 5 years ago (06/20/2016)
Hi Sia,Note that for points normalized by their norm, minimizing Euclidean distance is is equivalent to maximizing cross-correlation (aka "vector cosine"):\|x-y\|^2...
In Matlab or Octave, you can use freqz() to look at your FIR filter magnitude response.  You can set the lowpass cut-off frequency low enough so that its stop-band...
The signal-to-noise ratio in the surviving band should not change provided (1) your FFT length does not change its duration in seconds (giving same number of cycles...

Re: Newbie with questions about decimation for SDR

Reply posted 5 years ago (05/16/2016)
Also be aware of the sndfile-resample command-line program in the libresample package.  I use that all the time.  It is free open-source so you can study the software implementation. That...

Re: 90 Degrees Shift of Digital Signal

Reply posted 5 years ago (05/04/2016)
Here is a paper on designing a filter that will give phase quadrature over a wide frequency range:@conference{abel2010an, title = {An Infinite Impulse Response...

Re: 90 Degrees Shift of Digital Signal

Reply posted 5 years ago (04/28/2016)
By filtering a sinusoid to shift it 90 degrees, you are effectively already doing a Hilbert transform at that frequency. You are using a real pole to provide 45-degrees...

Re: Learning DSP with Python

Reply posted 6 years ago (04/04/2016)
I still use Octave and Matlab for high-level analysis scripting, but I contributed to a free online course that uses Python for its audio signal processing examples:https://www.coursera.org/course/audioI...

Re: s_to_z (Pupalaikis)

Reply posted 6 years ago (03/22/2016)
I have an old .m file for this kind of thing: deemphasis.m

Use this form to contact JOS

Before you can contact a member of the *Related Sites:

  • You must be logged in (register here)
  • You must confirm you email address