### Julius Orion Smith III (@JOS)

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

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

## Re: Narrow band filter at high rate

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

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

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...

## Re: Suggestions on book for "Fixed Point DSP with C/C++ "

I can recommend this: https://cnx.org/contents/feIAPEwX@1.1:vvZSF2za@2/S...

## Re: Matlab: frequency response of fir2 filter doesn't correspond to the used magnitude response

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

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

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...

## Re: Removing frequency components and reconstruct time domain signal.

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: Power spectral density

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...

...

## Re: Mathematical associative property of Convolution Reverb does not seem to hold. What am I doing wrong?

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

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...

## Re: Calculating the minimum phase of a given magnitude response

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...

## Re: How to work out the dynamic range of a filter empiracally

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?

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

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"

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.

For a more detailed discussion, you might tryhttps://www.dsprelated.com/freebooks/sasp/Hilbert_...

## Re: Frequency resolution

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) ?

Negating z is a pi-rotation, so you get a high-pass filter.

## Re: Generalized Cross Correlation method not producing desired results

I would also plot W, and limit its variation using the threshold I mentioned.

## Re: Generalized Cross Correlation method not producing desired results

Maybe change it to W = 1./(1.0e-4 + abs(Sxy));

## Re: Generalized Cross Correlation method not producing desired results

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...

## Re: create FIR filter from s-domain equation using Matlab firpm or other method

Check out cfirpm: https://www.mathworks.com/help/signal/ug/fir-filte...

## Re: What is the significance of whitening "signal + noise" at the receiver side in a wireless communication system

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...

## Re: What is the significance of whitening "signal + noise" at the receiver side in a wireless communication system

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

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

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

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...

## Re: How to retrieve information from a short sequence of the sum of two cosine waves with close frequencies

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

Yes, tikz is the gold standard.  I use xfig myself.

## Re: Spurious when Fs/F is not an integer

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?

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,...

## Re: Multi-loop, Multi-input Multi-output(MIMO) loop stability

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...

## Re: Multi-loop, Multi-input Multi-output(MIMO) loop stability

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

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

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...

## Re: Inverse Fourier transform and discrete time fourier series

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

I've never seen that - interesting - complexity nicely between TIIR and plain FIR

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

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

@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

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?

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

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...

## Re: how to apply Laplace filter value to FFT Nyquist bin?

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...

## Re: how to apply Laplace filter value to FFT Nyquist bin?

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

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

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

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...

## Re: Does decimation of a signal result in a shift in the baseline noise?

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...

## Re: Does decimation of a signal result in a shift in the baseline noise?

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

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

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

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

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)

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