Julius Orion Smith III (@JOS)
Fitting Filters to Measured Amplitude Response Data Using invfreqz in Matlab
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
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
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...
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...
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...
You should provide a small complete test example instead of the incomplete real thing.16002000 should probably be 1600 2000
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...
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.
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...
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...
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...
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...
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...
Apologies for likely duplication, but the optimal window for frequency-resolution in stationary noise is the rectangular window. See the classic papers by Rife...
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...
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...
The amplitudes of the harmonic-distortion components are smaller than the fundamental amplitude for a "weak nonlinearity", but not in general. Consider the squaring...
For frequency shifting, filter out the negative-frequency component and multiply by a similar "complex sinusoid" tuned to the "difference frequency", thereby summing...
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...
Another common practice is to do linear interpolation based on the low-order wavetable address bits. Below, for example, is the interpolating sinusoid oscillator...
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.
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
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...
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...
@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...
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...
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...
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...
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,...
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...
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...
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...
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...
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...
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...
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...
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
