Smoothing a staircase I have some data from a position encoder, so naturally i want to estimate its speed.  However, the data is very quantized, so it's difficult to smooth enough to differentiate easily:

Each step level is about 70-140 data points long on average, so my usual tricks of savitzky-golay filtering aren't up to it.  I've played with piecewise splines as well, but not had much luck; there's a lot of ringing.

Is there a straightforward way to draw the obvious line through the trace?

Looks like your data is virtually free of noise. That, combined with a very high sampling frequency would mean that at the jumps the data is exactly at the threshold between two quantized values. Set up nodes at the middle points of the vertical jumps and construct splines that connect the nodes. The easiest is to just draw straight lines between successive nodes, which gives a piece-wise constant differential. I wonder if that is good enough, or if you already tried that. If you need the velocity real-time, this approach is problematic because occasionally you might have to wait for a new node for some time.

You can further low-pass filter the interpolated data. If you use a filter with an impulse response that is nowhere negative, such as a Gaussian function, then there will be no overshoot.

With linear interpolation, everywhere between successive nodes, the speed will be simply the position difference between the nodes divided by the time difference between the nodes. You can run the smoothing filter on that piece-wise constant speed data and the result will be the same as if you'd run it on the linear position ramps and then differentiated (associative property of convolution, as both differentiation and filtering are convolution).


As I've said in the (copious) comments, the signal modeled using

$$
y[t]=A_1(\sin(\omega_1 t+\phi_1)+A_2 \sin(\omega_2 t+\phi_2)+....+A_p \sin(\omega_p t+\phi_p) + z(t)
$$

is stationary.

Your statement i am not interested estimation of parameters,just i am concerned about spectral resolution seems like an oxymoron.  I think what you mean is that you're not concerned with estimating $A_k$ and $\phi_k$, but you are interested in $\omega_k$ for the case where $\omega_k \approx \omega_{k+1}$.

Can you please refine your question, based on this information and the discussion in the comments above?


Polyphase filter bank channelizer FFT or IFFT I am referring to the Xilinx application note Polyphase Filter Bank Channelizer. The receiver portion should have FFT but the pdf mentions IFFT in both Tx and Rx.
The paper by F. J. Harris, C. Dick and M. Rice, "Digital receivers and transmitters using polyphase filter banks for wireless communications," in IEEE Transactions on Microwave Theory and Techniques, vol. 51, no. 4, pp. 1395-1412, April 2003 also mentions the same.
How did Xilinx get proper results by implementing IFFT in both Tx and Rx?

Sorry, I didn't read the paper so this is just a guess.
FFT and IFFT are almost identical algorithms so it really makes no sense to code up both in an FPGA. Just use
$$F^{-1} (x) = F(x')/N, \\ F(x) = N \cdot F^{-1}(x)'$$
where N is the length of the FFT and "'" to conjugate complex operator.


Step response overshoot in IIR filters Is there an analytic expression by which it is possible to calculate overshoot of step-response of second order IIR digital filter? I have tried to search through Mitra's signal processing book and through some other materials, but wasn't able to find it.

If we have a second-order discrete-time system with complex conjugate poles at $z_{\infty}=re^{\pm j\theta}$, $0<r<1$, $0<\theta<\pi$ (i.e. with poles inside the unit circle and no poles on the real axis), then we have a transfer function

$$H(z)=\frac{1}{1-2r\cos \theta z^{-1} + r^2z^{-2}}$$

For the moment we assume there are no zeros away from the origin of the $z$-plane. The step response of such a system is given by

$$a(n)=\frac{1}{1-2r\cos \theta + r^2}\left [1+
\frac{r^{n+2}\sin[\theta(n+1)] - r^{n+1}\sin[\theta(n+2)] }{\sin\theta} \right ]\sigma(n) =$$
$$=a_{\infty}[1+d(n)]\sigma(n)\tag{1}$$

where $\sigma(n)$ is the Heaviside step function, $a_{\infty}$ is the value of the step response for $n\rightarrow\infty$, and $d(n)$ corresponds to the ringing of the step response. Consequently, the maximum overshoot of $a(n)$ corresponds to the maximum of $d(n)$. You could try with Mathematica to see if there is an analytic solution for the (largest local) maximum of $d(n)$.

If the system has zeros away from the origin of the $z$-plane, then $H(z)$ has a numerator polynomal $b_0+b_1z^{-1}+b_2z^{-2}$:

$$H(z)=\frac{b_0+b_1z^{-1}+b_2z^{-2}}{1-2r\cos \theta z^{-1} + r^2z^{-2}}$$

In this case the formula for the step response given by (1) still holds with
$$a^{\prime}_{\infty}=a_{\infty}\sum_kb_k$$ and with

$$d^{\prime}(n)=\frac{(d*b)(n)}{\sum_kb_k}=\frac{b_0d(n)+b_1d(n-1)+b_2d(n-2)}{\sum_kb_k}$$

It may be even more difficult to find an analytic solution for the maximum of $d^{\prime}(n)$.

So unfortunately I cannot give a plug-and-play formula for the maximum overshoot, but my intent was to give you some analytical results to experiment with, either with Mathematica, or by using some very simple numerical search strategies.

If everything else fails (Mathematica etc.), I would simply compute $a_{\infty}$ by

$$a_{\infty}=\frac{b_0+b_1+b_2}{1+a_1+a_2}$$

and then actually compute the filter output (e.g. with Matlab's $\tt{filter}$ function), feeding the filter with a step, until you have found the first maximum of the output signal $a_{max}$. The overshoot is then given by $\frac{a_{max}}{a_{\infty}}$. But this may be obvious to you anyway.


Meaning of a null coefficient at 0 Hz I am analysing a wave with a timespan $\tau$ and after applying the fft I obtain the following plot of the coefficients

I assume that the peak at $t=1$ followed by smaller coefficients means that the wave has just one frequency, which is near $\frac{1}{\tau}$. Nevertheless, I would not expect a zero coefficient for 0Hz. What can be causing that zero coefficient? Is my interpretation wrong?.

A zero coefficient at DC simply means the mean of the waveform is 0; bin 0 of the DFT is identical to calculating the mean value scaled by $N$.
Consider the general expression for the DFT:
$$X(k) = \sum_{n=0}^{N-1}x(n)e^{-j k \omega_n n}$$
with $\omega_n = 2\pi/N$ and $k = 0,1 \ldots N-1$
(This is a correlation of the waveform x(n) with every integer multiple of the tone with frequency $\omega_n$)
When $k=0$ this reduces to:
$$X(0) = \sum_{n=0}^{N-1}x(n)$$
While the mean of $x(n)$  is
$$\bar x = \frac{1}{N}\sum_{n=0}^{N-1}x(n) = \frac{1}{N}X(0)$$
This is how the DFT works fundamentally; for each index k in the DFT, the signal $x(n)$ is frequency translated by $-k \omega_n$, and then we take the average of that result (without dividing by $N$, so a scaled average). If any frequency component was at $k \omega_n$ prior to the frequency translation, the resulting average would be large, indicating the presence of that frequency component.
Note for just one bin of the DFT how this is the basic structure of a Digital Downconverter as shown in the block diagram below for $x(n)$ real, where the low pass filters (LPF) are basic moving average filters (such as a CIC filter). The DFT computes one time sample, and replicates this for every single bin, and supports real or complex input signals where with a complex input there would be 2 additional multipliers and 2 adders for a full complex multiplication. So to put that in perspective, a 1023 point DFT is 1023 digital down-converters all running in parallel, which makes the FFT algorithm very impressive! Every computation of the DFT is one time sample of the result, so to process as a streaming receiver the input would be shifted one sample and then a new output computed (overlap-add and overlap-save methods improve the efficiency further for such streaming applications but shifting by one sample and computing each output makes one bin of the DFT and the digital-downconverter receiver shown below equivalent).
One Bin of the DFT as a DDC (for case of real input)


Is it okay to interpolate between overlapping areas obtained using a sliding window? I am using a 50x50 pixel sliding window which moves left to right every pixel so that it covers a whole binary image which is about 2000x2000 pixels. In the window, the ratio value of black to white pixels is calculated. I have then been uploading the values into matlab as xyz coordinates with x and y being the centre of each window square and z as the ratio value and interpolating between them to create a meshgrid colormap to help visualise the image.
This is outside the realms of my expertise but I was wondering is it okay to interpolate between these values in matlab as they are not really discrete points but overlapping areas. The colormaps come out looking really good and seem to be a very good representation of the original binary image but it's more the theory that I am having trouble with and I am not having much luck searching.
Is it okay to interpolate overlapping sliding window values in this way using matlab's standard linear interpolation mesh and does anyone have any links to helpful resources regarding this?
Any help would be greatly appreciated!

A first suggestion is to use an odd-sized sliding window (say $w=51$), so that you can really affect the black/white proportion to the center pixel. An even length can generate a half-pixel bias.
And you have a  window sliding by one pixel, therefore you have a black/white proportion for each pixel (except on the borders of the large image, of dimension $W\times W$). The resulting picture will be of size $(W-w+1)^2$, about $1950\times 1950$ here, which already seems to have a large definition for most displays. Hence, I cannot really understand the need to interpolate for visual purposes.
Nevertheless, two contiguous windows share $w\times (w-1)$ (sides) or $(w-1)\times (w-1)$ (diagonal) pixels in common, so more than 98% or 96 % respectively. Hence, the variations in the proportion should vary quite slowly. Then simple interpolation should work nicely. The problem is to give a meaning to the interpolated values. It should be close to having a window with fractional offsets over the big image, and counting partial pixels.


Isolation of a signal in a frequency band I have a continuous signal. I digitize it at 40kHz. A useful signal for me is in the range of 500-1500Hz and takes ~ 50-100Hz. It can be located in any part of this strip at any time.

If I apply a digital bandpass filter at 500-1500Hz, I will select the signal I want from the entire band, and along the way I will collect all the noise in the 500-1500Hz band. To reduce noise, I divided this band into 10 bands of 100 Hz (also with digital bandpass filters) and from them I choose the one with the largest signal. At the same time, I do not collect noise across all bands, that is, I have improved the signal-to-noise ratio.

But there are moments when the signal is on the border of my 10 filters, that is, part of it falls on one filter, part on another, and I lose signal-to-noise ratio.

An obvious option for me is to split the 500-1500Hz band into more bands, but the computational capabilities of my microcontroller are limited, and I cannot afford to increase the number of filters anymore.
Is there any good way to select a signal located in the 500-1500Hz band without collecting noise from the entire frequency band? I also wonder if there are 2 such signals in my 500-1500Hz frequency band, how can I do in this situation? And will it be possible to unite them? Perhaps the Wavelet transformation will help me, but in the Internet I did not find any information on how it can be implemented on a microcontroller.

The algorithm I referenced is premised on the signal being a pure real tone, in which case it will give an exact answer.
The presence of noise or other tones will distort the answer somewhat, but the formulas are actually quite robust.  As long as the tones are fairly far apart (at least several bins), the distortion due to other tones is small.
Aligning your DFT frame so the main tone is roughly a whole number of cycles and a half has two advantages.

The parameters are read from four values of the DFT instead of two, so inherently more noise resistant.

The second harmonic will fall near a whole number of cycles and thus not have much leakage.

If necessary, your signal can be filtered (as you have done) to reduce the magnitude of the harmonics, but you need to know the frequency response of the filter to adjust your amplitude calculation.
Real world signals tend not to be pure tones.  The smaller the segment, the more it will resemble a pure tone.  The formulas I referenced (and am the discoverer of) have the advantage of being exact, so a large N is not necessary for better results.  In a recent job, I was able to take 8 point DFTs on ~1 1/2 cycle frames to calculate the frequency, phase, and amplitude to about 4-5 significant figures.
In radar jobs in the past, and the non-radar one just mentioned, the variation in the signal meant the results you get are dependent on your frame selection and how you choose to average your results.  The algorithm can give you nearly instantaneous values along your signal duration.
The derivation of the frequency formula can be found here:

Two Bin Exact Frequency Formulas for a Pure Real Tone in a DFT

This is an extension to three bins, which can be used if your peak is near a bin (< 0.15 bin width) for more noise robust results.

Improved Three Bin Exact Frequency Formula for a Pure Real Tone in a DFT

In the pure noiseless single tone case, either formula will yield an exact answer.


What are some good signal and image processing blogs worth following? Are there some good sites or blogs where I can keep myself updated on the latest news and papers about image and signal processing research, or I should just check out "classical" providers like IEEE Transactions, Elsevier, etc?

There are many for different subjects -

Efg's algorithm collection : http://www.efg2.com/Lab/Library/ImageProcessing/index.html
DSP Forum : http://www.dsprelated.com/
Data compression - http://datacompression.info/
About rendering - http://www.realtimerendering.com/portal.html
For all research papers - http://ieeexplore.ieee.org/Xplore/guesthome.jsp
Resources on Mp3 and Audio - http://www.mp3-tech.org/programmer/docs/index.php
Steve on Image Processing - http://blogs.mathworks.com/steve/
Image Processing and Retrieval http://savvash.blogspot.com/
Accelerated Image Processing - http://visionexperts.blogspot.com/
The Digital Signal Processing Blog - http://centerk.net/dspblog/
Noise & Vibration Measurement Blog - http://blog.prosig.com/
Image Processing with Matlab, Open Blog - http://imageprocessingblog.com/


sampling rate for band pass signal The sampling frequency for band pass signal which is having frequency from 4 kHz to 6 kHz(rectangular shape from 4 to 6 kHz), can I prefer the sampling frequency of 12 kHz(2*the highest frequency present is message signal) or 6 kHz  (2*bandwidth of message signal(6-4)kHz), can u please suggest which one to prefer?

Sampling is frequency domain convolution by a frequency domain impulse train with the impulses spaced apart by the sampling frequency. If there is overlap between the images thus created, then there will be information loss. As you can see in the illustration below, all of the sampling frequencies $f_s$ = 12 kHz, $f_s$ = 6 kHz, and $f_s$ = 4 kHz will work fine. I tilted the band tops to indicate positive and negative frequencies. The first curve is the signal's spectrum before sampling. The unit is kHz.

Using a small sampling frequency has the benefit that any processing will have less samples that need to be computed per time unit. For example filters designed to a performance specification can be lower-order than with a higher sampling frequency.

Sometimes it is useful if the signal is oversampled though. Things like polynomial interpolation work better because of the more dense sampling. A Hilbert transformer would have more room for the transition band. If processing is non-linear, it may be possible to avoid aliasing of harmonics by filtering after the non-linear processing; with critical sampling the harmonics might alias and contaminate the signal.

Using a low sampling frequency (6 kHz or 4 kHz) that aliases the bands may require more work later if you need to have the band back at the original frequency.

So the optimal sampling frequency depends on what you are going to do with the signal.


Equalizing a temporal sound signal with a convolution This question is more or less the following of a previous question.

The aim is to equalize some sounds signals. By equalizing, I mean increasing or decreasing some frequency ranges, a bit like an analogic equalizer. The process I am applying now is :

1/ getting the temporal sound signal

2/ getting the "spectral correction" coefficients for each frequency sample

3/ converting the spectral correction to a "temporal correction" by using a FFT

4/ convoluting the temporal sound signal with the "temporal correction".

I have many questions with that process.

The first one is : what about the phasis ? The temporal sound frames are integers given by by wav file. So, if I write them with a complex number, the real part will be the value of the sample and the imaginary part will be 0.

But when i convert the spectral correction to a "temporal correction", some imaginary values may appears. And, after my convolution, I'll get a temporal signal with some imaginary parts wich will be different from 0. So, how, coming back to a temporal signal with a phasis = 0 ? Shall I use the modulus of the complex get by my convoluted signal ?

I also have a question about the "temporal correction". For testing, I'm using a "blank correction" : all the spectral coefficients = 1 so, my signal convoluted should be equal to the original signal. When I calculate the FFT of the spectral correction, I get a tab with those values {1, 0, 0, 0 (...)}.

So, my second question is : Does this result is logical or shall I get {1, 0, 0 (...), 1} because of the hermitian symmetry ?

And my third (and last) question is : What about the "padding" ? I'm following this website to compute the convolution but, I don't understand how many "0" I do have to add before or after the input signal. Let's say I have a 512 frames sound signal to convolute with a 511 frames "correction signal". Shall I add 511 "zero" before the signal ? or 255 before and 255 after ?

Thank you for all your answers.

Dr_Click

First, as long as your filter FFT has hermitian symmetry, the convolution with a real signal should also be real.  You may get some residual noise in the imaginary, but you can just pitch it.  If you are getting signal in the imaginary, your implementation is goofed.

Second, what you see is correct.  Because the time signal is real, the frequency signal has symmetry.  The opposite is also true.  Sample zero is the point about which the symmetry occurs.

Third, you have two lengths to take into account.  The sample window, N, and the filter kernel, M.  The sample window will determine how much real time lag you get through the operation, but is inconsequential in rendered applications.  The kernel length determines the number of “taps” in the filter bank.  The padding should be such that each signal is M+N-1 long when the operation is performed.

Fourth, this implementation is flawed because you are directly manipulating the signal in the frequency domain.  What you would want to do is something like this:  every N samples get M samples of the signal.  Perform an FFT.  Do whatever you need to to calculate the filter frequency response.  Set the phase to be linear and centered about the mid point.  Calculate the inverse FFT to get the filter kernel in the time domain.  Covolve the signals (w/ FFT convolution or directly).

Lastly, as I stated in response to your last question, this is not a good implementation of an analog prototype.  Analog filters have poles and should be implemented using biquad filters.  See RBJs biquad filter cookbook for a great source on implementing analog prototypes of common audio filters.  You can still do it the way I’ve described, but to relate it to the analog domain is a misnomer.


Fourier coefficients of discrete difference of a square wave I have a discrete square wave $f(t)$ where $t \in \mathbb{N}$, of amplitude $A$, period $T$ and duty cycle $1/T$
$$
f(t) = \left\{\begin{matrix}
A, & \mathrm{if}\;t=Tn\\
0, & \mathrm{if}\;t\neq Tn
\end{matrix}\right.
\;,\; n \in \mathbb{N}
$$
size = 23
T = 5
A = 10
t = np.arange(size)
f = np.zeros(size)
f[t%T==0] = A
fig, ax = plt.subplots(figsize=(12, 5))
ax.step(t, f)
ax.set(
    xticks=t,
    xlabel='time $t$',
    ylabel='amplitude $A$',
    title=f"Square wave of amplitude {A}, period {T} and duty cycle 1/{T}"
)
plt.show()

Given the Fourier series
$$
\mathcal{F}[g(t)] = a_0
+ \sum_{h=1}^{\infty} a_h \cos(2 \pi h \nu \cdot t)
- \sum_{h=1}^{\infty} b_h \sin(2 \pi h \nu \cdot t)
$$
for this square wave
$$
a_0 = \frac{A}{T}
$$
$$
a_h = \frac{2A}{\pi h} \sin \left( \frac{\pi h}{T} \right)
$$
$$
b_h = 0
$$
so
$$
\mathcal{F}[f(t)] = \frac{A}{T} + \sum_{h=1}^{\infty} \frac{2A}{\pi h} \sin\left(\frac{\pi h}{T}\right)
\cos\left(\frac{2 \pi h}{T} \cdot t\right)
$$
where the Fourier coefficients of harmonic $h$ is $\frac{2A}{\pi h} \sin\left(\frac{\pi h}{T}\right)$.
The discrete difference of $f(t)$
$$
f'(t) = f(t) - f(t-1)
$$
is the sum of two square waves $f(t)$ and $-f(t-1)$ with opposite amplitude and shifted by a sample unit
$$
f'(t) = \left\{\begin{matrix}
+A, & \mathrm{if}\;t=Tn\\
-A, & \mathrm{if}\;t=Tn+1\\
0, & \mathrm{otherwise}
\end{matrix}\right.
\;,\; n \in \mathbb{N}
$$
t1 = t[1:]
f1 = np.diff(f)

fig, ax = plt.subplots(figsize=(12, 5))
ax.step(t1, f1)
ax.set(
    xticks=t,
    xlabel='time $t$',
    ylabel='amplitude $A$',
    title=f"Discrete difference of a square wave of amplitude {A}, period {T} and duty cycle 1/{T}"
)
plt.show()

We can thus say (is it correct?) that the Fourier series of $f'(t)$ is the sum of the Fourier series $\mathcal{F}[f(t)]$ and $\mathcal{F}[-f(t-1)]$
$$
\mathcal{F}[f'(t)] =
\left[
\frac{A}{T} + \sum_{h=1}^{\infty} \frac{2A}{\pi h} \sin\left(\frac{\pi h}{T}\right)
\cos\left(\frac{2 \pi h}{T} \cdot t\right)
\right]
-
\left[
\frac{A}{T} + \sum_{h=1}^{\infty} \frac{2A}{\pi h} \sin\left(\frac{\pi h}{T}\right)
\cos\left(\frac{2 \pi h}{T} \cdot (t-1)\right)
\right]
$$
$$
\mathcal{F}[f'(t)] =
\sum_{h=1}^{\infty} \frac{2A}{\pi h} \sin\left( \frac{\pi h}{T}\right)
\left[
\cos\left( \frac{2 \pi h}{T} \cdot t\right) -
\cos\left( \frac{2 \pi h}{T} \cdot (t-1)\right)
\right]
$$
How can I get the Fourier coefficient of harmonic $h$ from $\mathcal{F}[f'(t)]$?
Is there a better way to express the Fourier series of $f'(t)$?

EDIT #1
As @MattL. noticed (see comments) what I'm dealing with is actually a discrete-time signal, so a better representation of the signal would be

and for the discrete difference will look like

It's quite straightforward to express the Fourier coefficients of a shifted version of a function in terms of the Fourier coefficients of the original function.
If the Fourier series of a $T$-periodic function $f(t)$ is given by
$$f(t)=\sum_{k=-\infty}^{\infty}c_ke^{j2\pi kt/T}\tag{1}$$
then it follows that the Fourier series of a shifted version must be
$$\begin{align}f(t-t_0)&=\sum_{k=-\infty}^{\infty}c_ke^{j2\pi k(t-t_0)/T}\\&=\sum_{k=-\infty}^{\infty}c_ke^{-j2\pi kt_0/T}e^{j2\pi kt/T}\\&=\sum_{k=-\infty}^{\infty}\tilde{c}_ke^{j2\pi kt/T}\tag{2}\end{align}$$
Hence, the Fourier coefficients of $f(t-t_0)$ are given by
$$\tilde{c}_k=c_ke^{-j2\pi kt_0/T}\tag{3}$$
And, by linearity, the Fourier series coefficients of $g(t)=f(t)-f(t-t_0)$ are
$$d_k=c_k-\tilde{c}_k=c_k\big(1-e^{-j2\pi kt_0/T}\big)\tag{4}$$

For periodic discrete-time sequences we can use the Discrete Fourier Tansform (DFT), which is just a discrete Fourier series representation:
$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{5}$$
with
$$X[k]=\sum_{k=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{6}$$
In a completely analogous manner to the continuous-time case, we can derive an expression for the DFT coefficients of the sequence $y[n]=x[n]-x[n-n_0]$, $n_0\in\mathbb{Z}$:
$$Y[k]=X[k]\big(1-e^{-j2\pi n_0k/N}\big)\tag{7}$$


Precise 5th and 7th harmonics of a sampled sine wave Does anyone know in decibels (to 1/100th of a dB) what the theoretical 3rd, 5th and 7th harmonic of a 0dB fs 24-bit (i.e. full-level; 0dB = -8,388,607 to 8,388,607) sampled sine wave without dither will be? Where the harmonics are a product of the sampling process (caused by the effect of the discrete sampled steps). Better still, how I calculate the theoretical values for the harmonic components.

The fundamental is around 1,000Hz (happy to use exactly 1kHz). The sampling frequency is 48kHz. Asynchronous sampling.

This answer discusses the harmonic spectra of the quantized sequence in five cases:

limit $f/f_s \to 0$,
synchronous sampling of a cosine with rational $f/f_s$,
synchronous sampling of a sinusoid of arbitrary phase with rational $f/f_s$,
asynchronous sampling of a sinusoid with rational $f/f_s$, and
asynchronous sampling of a sinusoid with rational $f/f_s$ in the limit $L/A \to
0$,

where $f$ is the frequency of the sinusoid, $f_s$ is the sampling frequency, $L$ is the fundamental period of the quantized sequence, and $A$ is the amplitude of the sinusoid. Decibels full-scale (dBFS) shall be referenced to the sinusoid amplitude $A$ at 0 dBFS.

Limit $f/f_s \to 0$

Let the zero-to-peak amplitude of the continuous-time prototypical sinusoid, which is being sampled and quantized, be a positive real number $A$. In the limit $f/f_s \to 0,$ where $f$ is the sinusoid frequency and $f_s$ is the sampling frequency, the amplitudes of the harmonics are found in the Fourier series representation of the continuous-time piece-wise constant quantized waveform. The terms of the Fourier series are directly the harmonics of the fundamental frequency. We can construct the waveform as a sum of components like the red curve here, illustrated for sinusoid amplitude $A = 7$:

Figure 1. prototypical sinusoid, of amplitude $A=7$, and its quantization in blue, with a component highlighted in red.

We are working with cosine instead of a sine because the math is nicer this way. Because of the symmetries, each component has only odd harmonics and only cosine terms in its Fourier series. The non-zero coefficients of the Fourier series of a component of integer amplitude $k \in 1\ldots \operatorname{round}(A)$ are given by:

$$a_n = \frac{4}{\pi}\int_{x_0}^{x_1}\cos(nx)\,k\,dx = \frac{4k}{\pi n}\big(\sin(x_1n) - \sin(x_0n)\big),\quad n\text{ odd},\tag{1}$$

with the component having the said amplitude over the interval $x_0 < x < x_1$ in the first quarter-period of the cosine. I did not bother to explicitly write the contributions from the symmetrically arranged non-zero pieces of the component in the other quarters, because they will contribute identically to the odd harmonic cosine terms in the Fourier series. Instead I simply included the implicit factor 4 in the equation.

The Fourier series of the full piece-wise constant waveform is a sum of the Fourier series of the components. The boundaries of the pieces that we need to include in the sum are:

$$\begin{gather}0 < x < \operatorname{acos}\left(\frac{\operatorname{round}(A) - 0.5}{A}\right),\quad \text{if }k = \operatorname{round}(A),\\
\operatorname{acos}\left(\frac{k + 0.5}{A}\right) < x < \operatorname{acos}\left(\frac{k - 0.5}{A}\right),\quad \text{if }k \in 1\ldots \operatorname{round}(A)-1,\end{gather}\tag{2}$$

where $\operatorname{acos}\left(\frac{s}{A}\right)$ comes from solving $A\cos(x) = s$ in the first quarter-period. The non-zero coefficients of the Fourier series of the full piece-wise constant waveform are then given by:

$$\begin{gather}a_n = \frac{4\operatorname{round}(A)}{\pi n}\,\sin\bigg(\operatorname{acos}\left(\tfrac{\operatorname{round}(A) - 0.5}{A}\right)n\bigg) +\\\frac{4}{\pi n}\sum_{k=1}^{\operatorname{round}(A)-1} \Bigg(k\sin\bigg(\operatorname{acos}\left(\tfrac{k - 0.5}{A}\right)n\bigg) - k\sin\bigg(\operatorname{acos}\left(\tfrac{k + 0.5}{A}\right)n\bigg)\Bigg),\quad n\text{ odd}.\end{gather}\tag{3}$$

The rounding convention does not matter as long as it is to a nearest integer.

8-bit quantization

It is not out of reach to compute Eq. 3 for $A=2^{23} - 1$ and $n \in {5, 7},$ as parameterized in the question, using something like Python's mpmath. Let's try that but with a 8-bit sinusoid first, because we can then evaluate $a_n$ at multiple $n$ in a reasonable time and perhaps learn something from the results. In Python:

import numpy as np
import mpmath as mp
def a_n(A, n):
    if int(n) & 1:
        return 4*mp.floor(A + 0.5)/(mp.pi*n)*mp.sin(mp.acos((mp.floor(A + 0.5)-0.5)/A)*n) + 4/(mp.pi*n)*mp.nsum(lambda k: k*( mp.sin(mp.acos((k-0.5)/A)*n) - mp.sin(mp.acos((k+0.5)/A)*n) ), [1, mp.floor(A + 0.5)-1])
    else:
        return 0

A = 2**7-1  # amplitude (positive real)

for n in np.arange(1, 10**4, 2):
    print(str(n)+","+str(20*mp.log10(mp.fabs(a_n(A, n))/A)))  # n, a_n (dBFS)

Figure 2. Amplitudes of harmonics of the quantized prototypical sinusoid of amplitude $A=2^{7}-1$. Only odd harmonics are present.

It took some minutes on my PC to get the results (Fig. 2). The graph looks jagged. I think this is because the waveform is similar to a sawtooth riding on a sinusoid, and near the zero crossing the sawtooth frequency stays approximately the same over a large portion of the waveform, being $2\pi A \approx$ 797.964534 times the fundamental frequency at the zero crossing as an asymptotic approximation in the limit $A\to\infty$ of the frequency $\pi/\operatorname{asin}\big((k + 0.5)/A\big)$ calculated based on the distance between the discontinuities of the saw wave just around the zero crossing. Indeed, the largest amplitude is found nearby at $|a_{787}| =$ -68.30 dBFS. This is much higher than the amplitude of a "random" early harmonic. Such peaks are repeated later at approximate harmonics of the sawtooth frequency but less clearly and with a roll off arising from the factor $1/n$ in Eq. 3, making the asymptotic decay -20 dB/decade (~-6 dB/octave) as $n\to\infty$. The amplitude of the fundamental frequency is 0.00066 dBFS, I think because the rounding at the top of the sinusoid is "pulling it up" from the 0 dBFS amplitude of the prototypical cosine.

24-bit quantization

The following continuation of the earlier Python script computes in high precision your desired numbers for $A=2^{23} - 1$, in the limiting case $f/f_s \to 0$:

A = 2**23-1 # amplitude (positive real)

n = 5  # number of the harmonic (odd positive integer)
mp.mp.prec = 70  # Precision 1
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dBFS
mp.mp.prec = 140  # Precision 2
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dBFS

n = 7  # number of the harmonic (odd positive integer)
mp.mp.prec = 70  # Precision 1
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dB
mp.mp.prec = 140  # Precision 2
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dB

mp.mp.prec = 53  # Restore default precision

I'm only reporting as many digits as are agreed about by repeated computation using two different mpmath precision settings, including also the digits that got changed by $\pm1$ due to larger changes in later digits. After a total of a few hours of computation on my PC, the result is that with $A=2^{23} - 1$ and $f/f_s \to 0,$ the amplitude of the 5th harmonic is -226.91150085 dBFS and the amplitude of the 7th harmonic is -226.9115030 dBFS. I also computed the amplitude of the fundamental frequency, 3.9195785E-11 dBFS. I also computed in a full day the amplitudes of a number of higher harmonics, in Python:

A = 2**23-1 # amplitude (positive real)

for m in range(1, 20):
  n = mp.mpf(10)**m+1 # number of the harmonic
  mp.mp.prec = 70  # precision
  a_n_low_prec = a_n(A, n)
  mp.mp.prec = 90  # precision
  a_n_high_prec = a_n(A, n)
  print(str(n)+","+str(20*mp.log10(mp.fabs(a_n_low_prec)/A))+","+str(20*mp.log10(mp.fabs(a_n_high_prec)/A)))  # n, a_n

Figure 3. Amplitudes of a selection of odd harmonics of the quantized prototypical sinusoid, with $A=2^{23}-1$.

The corner frequency before the roll-off is again approximately the sawtooth frequency at the zero crossing of the quantized prototypical sinusoid, $2\pi A = 2\pi(2^{23}-1) \approx$ 5.270717225E7 times the fundamental frequency $f$.

Synchronous sampling of a cosine with rational $f/f_s$

Let the prototypical sinusoid, which is being sampled and quantized, be a cosine without a phase shift. That is to say, one of the samples is at the peak of the sinusoid. Let $f/f_s$ be equal to a rational number $c/d$, with $c$ and $d$ integer. The fundamental period $L$ of the periodic sequence of quantized samples is:

$$L = \frac{d}{\operatorname{gcd}(c, d)},\tag{4}$$

where gcd denotes the greatest common divisor. A real discrete Fourier transform (real DFT) of length $L$ has bin frequencies that are the harmonic frequencies of the period-$L$ frequency. The period-$L$ frequency is equivalent to the frequency of the prototypical sinusoid only if the numerator of the irreducible fraction representation of $f/f_s$ equals 1. Otherwise the real DFT bin frequencies contain also other frequencies than those equivalent to in-band harmonics of the sinusoid.

What we have learned from the results so far is that the amplitudes of the harmonics of the quantized prototypical sinusoid do not start to roll off until after the approximately $2\pi A$th harmonic. This means that with rational $f/f_s$, the amplitude of an in-band harmonic frequency of the quantized sequence is typically almost fully determined by aliased frequencies. As an example of such aliasing, with $f_s =$ 48 kHz and $f =$ 1000 Hz, the 5th harmonic gets summed with the aliases of the 43th, 53th, 91th, 101th, 139th, 149th, etc. harmonic:

$$\begin{eqnarray}
&& \ldots\\
&=& 48000\text{ Hz}\times -3 + 1000\text{ Hz}\times 149\\
&=& 48000\text{ Hz}\times -2 + 1000\text{ Hz}\times 101\\
&=& 48000\text{ Hz}\times -1 + 1000\text{ Hz}\times 53\\
&=& 1000\text{ Hz}\times 5\\
&=& 48000\text{ Hz}\times 1 - 1000\text{ Hz}\times 43\\
&=& 48000\text{ Hz}\times 2 - 1000\text{ Hz}\times 91\\
&=& 48000\text{ Hz}\times 3 - 1000\text{ Hz}\times 139\\
&=& \ldots\end{eqnarray}\tag{5}$$

The effect of summation of the harmonics aliased to a frequency depends also on the phase of the sinusoid being quantized, because the phases of the aliased harmonics control whether there is constructive or destructive interference. Aliasing does not change the phase of a harmonic. A phase shift of the quantized prototypical sinusoid is equivalent to shifting it in time, which accounts to a phase shift of each unaliased harmonic proportional to its number. In this answer we will only be sampling and quantizing a prototypical cosine without a phase shift, so to compute a single bin of the real DFT of the discrete-time periodic sequence of samples, it would suffices to sum those $a$ coefficients of the Fourier series of the quantized prototypical sinusoid that alias to the equivalent frequency of the bin. We'll use the normalization convention that $|\operatorname{bin}_1| = 1$ for the real DFT of a sinusoid of amplitude $1$ and of angular frequency $2\pi/L$, where $L$ is the length of the real DFT in time domain. If the numerator of $f/f_s$ expressed as an irreducible fraction equals 1, we have:

$$\operatorname{bin}_n = a_n + \sum_{k=1}^{\infty}\left(a_{kL-n} + a_{kL+n}\right),\tag{6}$$

where $a_n$ is given by Eq. 3 for odd $n$ and $a_n$ is zero for even $n$. With $f = $ 1000 Hz, $f_s = $ 48 kHz, on my PC it would take perhaps thousands of years to compute Eq. 6 to reasonable accuracy for $A=2^{23}-1$. I also have difficulties evaluating the series even with $A=2^3-1$, with the sequence of partial sums oscillating. I think this difficulty arises from that with $A = 2^3 - 1$ and $f/f_s = 1/48$, the sequence before quantization contains values $-3.5$ and $3.5$ exactly at a discontinuity in the quantized prototypical sinusoid. In any case, with rational $f/f_s$, the real-DFT-based approach to determine the amplitudes of the harmonics is much more practical.

24-bit quantization

The following Python script computes using the real-DFT-based method the amplitudes of the harmonics of a synchronously sampled cosine for $A = 2^{23}-1$, $f =$ 1000 Hz, and $f_s =$ 48 kHz:

import numpy as np
from sympy import Rational
def bins(A, c, d):
    L = Rational(d, np.gcd(c, d))  # Fundamental period (samples)
    waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d)))
    return np.fft.rfft(waveform)/(A*L/2)

A = 2**23-1  # amplitude (positive real)
c = 1000  # numerator of f/f_s (positive integer)
d = 48000  # denominator of f/f_s (positive integer)
[20*np.log10(float(abs(x))) for x in bins(A, c, d)]  # real DFT magnitudes (dBFS)

While not required by the script, in this case $c/d = 1/48,$ so each real DFT bin corresponds to a harmonic of the prototypical sinusoid. The fundamental frequency's amplitude measures as -1.75E-7 dBFS, the 5th harmonic as -160.90 dBFS and the 7th as -160.75 dBFS. Unfortunately I could not find a suitable multiple precision Fast Fourier Transform (FFT) library, but I think those three numbers came out in sufficient precision using NumPy's rfft. The precision is not sufficient to calculate the amplitude of say the 3rd harmonic seemingly at about -333 dBFS (Fig. 4).

Figure 4. Magnitudes of the real DFT bins of the quantized sequence of a synchronously sampled cosine with $A=2^{23}-1$, $f =$ 1000 Hz, and $f_s =$ 48 kHz. These results are only applicable when one of the samples is at the peak of the sinusoid, so not applicable to arbitrary phase shifts or asynchronous sampling.

With this $f/f_s$ the fundamental period is $L = 48$. The computed amplitudes of the harmonics are also correct for phase shifts of the prototypical cosine that are multiples of $2\pi/L$, because the cosine is sampled also at those phases, or such a multiple plus $\pi$, because of the symmetry properties of cosine. In the presented case the compatible phase shifts of the cosine are multiples of $7.5°$.

We have encountered here a troublesome combination of $L = 48$ and $A = 2^{23}-1$. We have samples $A\cos(8\times2\pi/L) =$ $A\cos(40\times2\pi/L) =$ $A/2$ and $A\cos(16\times2\pi/L) =$ $A\cos(32\times2\pi/L) =$ $-A/2$, and as $A$ happens to be an odd integer, those samples will always have a fractional part of $1/2$. The above reported amplitudes of the harmonics may thus be dependent on the rounding mode, or, the computation of the phase and/or evaluation of the cosine may have given a small numerical error which has determined the direction of rounding. This may even break the symmetry properties of the sequence, affecting which amplitudes of the harmonics of the period-$L$ frequency are zero.

Method validation

We can validate the two approaches against each other by choosing a pair of $A$ and $f/f_s$ that does not result in sampling of the quantized prototypical cosine at any of its discontinuities. Firstly, this removes having to choose arbitrarily which way a number exactly half-way between integers should be rounded, and secondly, should help in convergence of the infinite series in Eq. 6. Let's try $A = 8$ and $f/f_s = 1/48$. First, the real-DFT-based method, in Python:

import numpy as np
from sympy import Rational
def bins(A, c, d):
    L = Rational(d, np.gcd(c, d))  # Fundamental period (samples)
    waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d)))
    return np.fft.rfft(waveform)/(A*L/2)

A = 8  # amplitude (positive reals)
c = 1000  # numerator of f/f_s (positive integer)
d = 48000  # denominator of f/f_s (positive integer)
[20*np.log10(float(abs(x))) for x in bins(A, c, d)]  # real DFT magnitudes (dBFS)

This results in the following amplitudes of the harmonics, printed in precision beyond their accuracy, in dBFS:

-inf, 0.005622747208892056, -inf, -331.19944825653926, -inf, -50.41376796795221, -inf, -35.41599672115829, -inf, -327.2878768807577, -inf, -38.14783193548751, -inf, -45.42857299685606, -inf, -324.8188873214554, -inf, -40.41625528549065, -inf, -36.12873038111221, -inf, -337.55682735433885, -inf, -33.04816436790989, -inf

Then, aliasing the harmonics of the quantized prototypical cosine by evaluating Eq. 6, with the series truncated to the first approximately $10^4\pi A$ terms, continuing with the imports and constants of the previous Python script:

import mpmath as mp
def a_n(A, n):
    if int(n) & 1:
        return 4*mp.floor(A + 0.5)/(mp.pi*n)*mp.sin(mp.acos((mp.floor(A + 0.5)-0.5)/A)*n) + 4/(mp.pi*n)*mp.nsum(lambda k: k*( mp.sin(mp.acos((k-0.5)/A)*n) - mp.sin(mp.acos((k+0.5)/A)*n) ), [1, mp.floor(A + 0.5)-1])
    else:
        return 0

mp.mp.prec = 100
L = Rational(d, np.gcd(c, d))  # Fundamental period (samples)
if c == np.gcd(c, d):
    for n in range(L/2):
        bin = a_n(A, n) + mp.nsum(lambda k: a_n(A, k*L - n) + a_n(A, k*L + n), [1, int(10000*mp.pi*A)])
        float(20*mp.log10(mp.fabs(bin)/A)) # harmonic amplitude (dBFS)
else:
    print("Error: fundamental frequency must be in bin_1")

mp.mp.prec = 53

which results in the following amplitudes, printed in precision beyond their accuracy, in dBFS:

-inf, 0.005622452922818016, -inf, -165.35811090601203, -inf, -50.41370472648242, -inf, -35.41599394709451, -inf, -155.34167605811237, -inf, -38.14782298030069, -inf, -45.42859376902997, -inf, -164.6586753088176, -inf, -40.41623886311205, -inf, -36.12872396082294, -inf, -149.78722963456818, -inf, -33.04815743325275

The results seem to agree within their numerical accuracy (Fig 6.).

Figure 6. Comparison of amplitudes, in dBFS, of harmonics of a quantized cosine sampled synchronously, with $A = 8$ and $f/f_s = 1/48$, calculated by the real-DFT-based method (blue squares) and by summing coefficients $a_n$ of harmonics that alias to the same frequency (orange crosses).

The 3rd harmonic amplitude calculated by truncating the sum in Eq. 6 converges very slowly as function of the number of terms (Fig. 6). This is probably the case also with the other low-amplitude harmonics. It does not seem feasible to reach the accuracy of the real-DFT-based method using the harmonic aliasing method.

Figure 7. $\operatorname{bin}_3$ as function of the number of terms in truncation of the sum in Eq. 6.

Synchronous sampling of a sinusoid of arbitrary phase with rational $f/f_s$

24-bit quantization

The real-DFT-based method is simple to modify for arbitrary phase shifts of the prototypical cosine. For $A = 2^{23}-1$, $f =$ 1000 Hz, $f_s =$ 48 kHz, and a phase shift $\phi = 0.123$ in radians, the amplitudes of the harmonics of the fundamental period can be computed by the following Python script:

import numpy as np
from sympy import Rational
def bins(A, c, d, phi):
    L = Rational(d, np.gcd(c, d))  # Fundamental period (samples)
    waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d) + phi))
    return np.fft.rfft(waveform)/(A*L/2)

A = 2**23-1  # amplitude (positive real)
c = 1000  # numerator of f/f_s (positive integer)
d = 48000  # denominator of f/f_s (positive integer)
phi = 0.123  # phase shift of the cosine (radians)
[20*np.log10(float(abs(x))) for x in bins(A, c, d, phi)]  # real DFT magnitudes (dBFS)

Figure 8. Magnitudes of the real DFT bins of the quantized sequence of a synchronously sampled cosine with phase shift $\phi = 0.123$, $A=2^{23}−1$, $f=$ 1000 Hz, and $f_s$= 48 kHz.

Asynchronous sampling of a sinusoid with rational $f/f_s$

Asynchronous sampling with rational $f/f_s$ can be modeled by assuming that $f/f_s$ is very close to a rational number. Then, locally the sequence of samples appears periodic with a computable harmonic spectrum. After sufficient time has elapsed, the phase of the prototypical sinusoid will have drifted to a different value and with a different local harmonic spectrum. For a random such sufficient time, the drift is a random variable between $0$ and $2\pi$. A reasonable probabilistic analog of the amplitude of a real DFT bin of the periodic sequence is the expected value of the square of its absolute value, $\operatorname{E}[|\operatorname{bin}_n|^2]$, which can be expressed in the same dBFS scale. We can estimate its value using uniform statistical sampling over phase shifts $\phi$ in a non-redundant range $0 < \phi < 2\pi/L$, where $L$ is the fundamental period of the locally periodic sequence. The end points of the range shall be deliberately left out so that we don't have the previously encountered problem with rounding modes and such.

24-bit quantization

$\operatorname{E}[|\operatorname{bin}_n|^2]$ is estimated for $A = 2^{23}-1$, $f =$ 1000 Hz, and $f_s =$ 48 kHz, by the following Python script, which also estimates $\max|\operatorname{bin}_n|$ over random phase shifts, representing the worst-case non-fundamental harmonic amplitudes in the local spectra:

import numpy as np
from sympy import Rational
def E_sq_abs_bins_and_E_max_sq_abs_bins(A, c, d, N):
    L = Rational(d, np.gcd(c, d))  # Fundamental period (samples)
    num_bins = int((L + 2)/2)
    sum_sq_bins = np.zeros(num_bins)
    max_sq_bins = np.zeros(num_bins)
    for n in range(N):
        phi = float((0.5 + n/N)*2*np.pi/L)
        waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d) + phi))
        sq_bins = np.square(np.absolute(np.fft.rfft(waveform)/(A*L/2)))
        sum_sq_bins = sum_sq_bins + sq_bins
        max_sq_bins = np.maximum(max_sq_bins, sq_bins)
    return [sum_sq_bins/N, max_sq_bins]

A = 2**23-1  # amplitude (positive real)
c = 1000  # numerator of f/f_s (positive integer)
d = 48000  # denominator of f/f_s (positive integer)
N = 2**18 # size of statistical sample (positive integer)
results = E_sq_abs_bins_and_E_max_sq_abs_bins(A, c, d, N)
[10*np.log10(float(x)) for x in results[0]]  # est. expected value of real DFT squared magnitudes (dBFS)
[10*np.log10(float(x)) for x in results[1]]  # est. max real DFT squared magnitudes (dBFS)

The results convergence quite slowly as the size of the statistical sample N is increased, and the accuracy requirements stated in the question are almost but not quite fulfilled by a day's calculation on my PC for $\operatorname{E}[|\operatorname{bin}_n|^2]$ and not at all for $\max|\operatorname{bin}_n|$. In the results obtained (Fig. 9) $\operatorname{E}[|\operatorname{bin}_n|^2]$ is about -157.1 dB and $\max|\operatorname{bin}_n|$ about -147.2 to -146.0 dBFS for non-fundamental harmonics. For the fundamental frequency, $\operatorname{E}[|\operatorname{bin}_1|^2]$ is in the order of -0.95E-11 dBFS, and $\max|\operatorname{bin}_1|$ is in the order of 3.4E-7 dBFS or higher.

Figure 9. Estimated $\operatorname{E}[|\operatorname{bin}_n|^2]$ (blue squares) and $\max|\operatorname{bin}_n|$ (red lines) for local $\operatorname{bin}_n$ in slow-drifting asynchronous sampling, for $A = 2^{23}-1$, $f =$ 1000 Hz, and $f_s =$ 48 kHz, estimated by a uniform statistical sample of phase shifts with a sample size of $2^{18}$.

The estimates show that asynchronous sampling gives a smoothing effect on the statistics of the local amplitudes of the harmonics, compared to synchronous sampling with a fixed alignment of the prototypical cosine with respect to the sampling grid. However, with asynchronous sampling, sometimes the worst case local harmonic peaks occur, whereas with synchronous sampling the harmonic amplitudes are well controlled and the worst cases are avoided by a suitable choice of phase shift.

Asynchronous sampling of a sinusoid with rational $f/f_s$ in the limit $L/A\to 0$

In this section of the answer, we will consider only the question's case that $L$ is even and $f/f_s = 1/L$, that is to say that $\operatorname{bin}_1$ corresponds to the frequency of the prototypical sinusoid.

If the length of the fundamental period of the locally periodic sequence of samples $L$ calculated by Eq. 4 is very small compared to the amplitude $A$ of the prototypical sinusoid, then for different phase shifts of the sinusoid, the quantization error can perhaps be considered a sequence of independent pseudorandom numbers distributed uniformly from -0.5 to 0.5, with an antiperiodicity constraint if $L$ is even (and some other constraint if $L$ is odd). Independence of the quantization errors would require that at a number of certain small time lags the autocorrelation of the quantized prototypical sinusoid approaches zero, but this is not something we know to be a fact. If it was true, then in the limit $L/A \to 0$ any sequence consisting of ±0.5's could be found as a quantization error sequence at some phase shift $\phi$ of the prototypical cosine, and within those sequences there would be sequences that maximize $|\operatorname{bin}_n|$ of the real DFT of the quantization error sequence. Due to orthogonality of the real DFT basis vectors, if $n$ does not correspond to the fundamental frequency of the sinusoid, then $\operatorname{bin}_n$ of the quantized sequence equals $\operatorname{bin}_n$ of the quantization error. If $n$ corresponds to the fundamental frequency, then $\operatorname{bin}_n$ of the quantized sequence equals $\operatorname{bin}_n$ of the quantization error plus $Ae^{i\phi}$.

Even if we do not assume independence of quantization errors, we can use the line of thinking of the previous paragraph to find an upper bound of $\max|\operatorname{bin}_n|$, the maximum absolute bin value for a harmonic of the prototypical sinusoid. A worst-case sequence of ±0.5's that maximizes $|\operatorname{bin}_n|$ could be found by an exhaustive search in the binary sequences, but can also obtained directly. A worst-case sequence is any sequence $\operatorname{seq}_k = 0.5\operatorname{sgn}\big(\cos(2\pi nk/L + \chi)\big)$, where $\operatorname{sgn}$ is the sign function, and with a phase shift $\chi$ that does not result in a zero crossing of the expression's sinusoid being at any sample point (the integers). One worst-case sequence is found in the limit of approaching from one side the transition between two worst-case sequences:

$$\begin{eqnarray}\operatorname{seq}_k &=& 0.5\times\begin{cases}\operatorname{sgn}\big(\sin(2\pi nk/L)\big),&\text{if }\sin(2\pi nk/L) \ne 0,\\
\cos(2\pi nk/L),&\text{if }\sin(2\pi nk/L) = 0,\end{cases}\\
&=& 0.5\times\begin{cases}\operatorname{sgn}\big(\sin(2\pi nk/L)\big),&\text{if }2kn \text{ mod } L \ne 0,\\
\cos(2\pi nk/L),&\text{if }2kn \text{ mod } L = 0,\end{cases}\end{eqnarray}\tag{7}$$

where $\text{ mod }$ denotes remainder in integer division in a numerically safe expression. An upper bound of $\max|\operatorname{bin}_n|$ is then found by real DFT of the worst-case sequence.

24-bit quantization

For odd $n$ the sequence given by Eq. 7 satisfies the antiperiodicity constraint for even $L$, and we can calculate an upper bound of $\max|\operatorname{bin}_n|$ of the maxima of local harmonic amplitudes $|\operatorname{bin}_n|$ for $A = 2^{23}-1$, $f =$ 1000 Hz and $f_s =$ 48 kHz, in Python:

import numpy as np
from sympy import Rational
def bins_max_err_seq(A, c, d, n):
    L = int(Rational(d, np.gcd(c, d)))  # Fundamental period (samples)
    waveform = 0.5*np.sign(np.sin(np.dot(range(L), 2*np.pi*n/L)))
    for k in range(L):
        if 2*k*n % L == 0:
            waveform[k] = 0.5*np.cos(k*2*np.pi*n/L)
    return np.fft.rfft(waveform)/(A*L/2)

A = 2**23-1  # amplitude (positive real)
c = 1000  # numerator of f/f_s (positive integer)
d = 48000  # denominator of f/f_s (positive integer)
if c == np.gcd(c, d) and int(Rational(d, np.gcd(c, d))) % 1 == 0:
    for n in np.arange(1, int(d/c/2 + 1), 2):
        print(str(n) + "," + str(20*np.log10(float(abs(bins_max_err_seq(A, c, d, n)[n])))))  # real DFT magnitude of bin n for the maximum quantization error sequence (dBFS)
else:
    print("Error: Fundamental frequency must be in bin 1 and L must be even\n")

For the odd harmonics, for $A = 2^{23}-1$, $f =$ 1000 Hz and $f_s =$ 48 kHz, we get locally $\max|\operatorname{bin}_n| \le $ -142.34 dBFS when the number of the harmonic $n$ is divisible by 3, and $\max|\operatorname{bin}_n| \le $ -142.39 dBFS otherwise which applies to the 5th and to the 7th harmonic. This seems compatible with Fig. 9, with an about 4 dB higher upper bound of the maximum local amplitudes of harmonics compared to the maximum local harmonic amplitude found by numerical search in actual quantization error sequences. A willing eye can also see an elevation in $\max|\operatorname{bin}_n|$ of Fig. 9 when $n$ is divisible by 3.

This answer has reached the maximum allowed length of 30000 characters.


Complex signal low pass filter I'm trying to implement a low pass filter in python for a complex signal but the output doesn't look right. I've created a simple example below where I've mixed a 15Hz complex sine wave and a 30Hz complex sine wave so that I get a signal with components (30+15=45Hz and 30-15=15Hz).
What I'm trying to do is filter this signal with a low pass filter to just give the resultant 15Hz component. When I look at the plot of the filtered signal though the imaginary component seems to be almost identical (same phase) as the real signal. It doesn't look right.
Example of code to reproduce then plot below:
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, lfilter
import scipy.signal as sig

#Sample rate.
rate = 5000

# Create a 15 Hz comlex signal.
t = np.arange(0,5000)/rate
f1 = 15
msg =  np.exp(1j* (2*np.pi * f1 * t))

# Create a 30 Hz complex signal.
f2 = 30
bb_lo = np.exp(1j * (2*np.pi * f2 * t))

# Create a low pass filter to remove the upper frequency.
nyq = 0.5 * rate
cut = 30
low = cut / nyq
b, a = sig.butter(5, low, btype='low')

# Mix the low frequency and high frequency signal together.
# Will result in a signal with freq components of f1+f2 and f2 - f1.
# Then apply the filter to real and imaginary parts.
signalFilt = lfilter(b,a , bb_lo.real*msg.real) + 1j*lfilter(b,a , bb_lo.imag*msg.imag)

plt.plot(signalFilt.real)
plt.plot(signalFilt.imag)
plt.grid()
plt.show()

Frequency filter looks right at approx 15Hz but real/imaginary in plot I don't think should be the same. What am I doing wrong here?
Thanks

You're not actually mixing two complex exponentials, because otherwise you'd get a complex exponential with the sum of the frequencies of the individual complex exponentials. That was already mentioned in a comment by Marcus Müller.
What you are doing is compute a signal with real part $\cos(\omega_0t)\cos(\omega_1t)$ and with imaginary part $\sin(\omega_0t)\sin(\omega_1t)$. If you lowpass filter that signal, only the components with the difference of the two frequencies remain, but these components are the same for the real and for the imaginary part, because
$$\cos(\omega_0t)\cos(\omega_1t)=\frac12\big[\cos[(\omega_0-\omega_1)t]+\cos[(\omega_0+\omega_1)t]\big]\tag{1}$$
and
$$\sin(\omega_0t)\sin(\omega_1t)=\frac12\big[\cos[(\omega_0-\omega_1)t]-\cos[(\omega_0+\omega_1)t]\big]\tag{2}$$
Clearly, the components with the difference frequency are the same.
Also, if you want to filter out a component with frequency $\omega_x$ (in your case the component with the sum frequency), it is not wise to use a filter with cut-off frequency $\omega_c=\omega_x$ because the attenuation at $\omega_c$ is usually quite modest, in your case $3$dB. Choose a frequency that is much closer to the lower frequency component to get a decent attenuation of the higher frequency component.


What's wrong with this implementation of my IIR filter? I posted this originally in the electronics stack exchange, as it's for an electronic purpose (digitally controlled boost) but it was pointed out to me that this would be the more appropriate place - so here goes!

As mentioned, I recently started developing a digital regulator for a DCM boost and I've implemented the following discrete transfer function in C:

$$ G_{Z} = \frac{1.39 z^{-1} - 1.332z^{-2} }{2 - 2.002{z^{-1}} + 0.002031z^{-2}} $$

and, in terms of a difference equation:

$$ y[n] = 0.695\cdot y[n-1] - 0.666 \cdot y[n-2] + 1.001\cdot x[n-1] - 0.0010155 \cdot x[n-2] $$

My code is as follows:

// Digi Reg Coeffs
const float a1 = 1.001;
const float a2 = 0.0010155;

const float b1 = 0.6975;
const float b2 = 0.666;

typedef struct {
  float y0; // new reg value
  float y1; // last reg value
  float y2; // last^2 reg value

  float x0; // latest sample
  float x1; // last sample
  float x2; // last^2 sample
} voltage_control;

uint16_t DigitalRegulator(float error) {
  float duty = 0;

  voltage_control.x2 = voltage_control.x1;
  voltage_control.x1 = voltage_control.x0;
  voltage_control.x0 = error; /// this is the new sample

  voltage_control.y2 = voltage_control.y1;
  voltage_control.y1 = voltage_control.y0;

  // calc new y0
  voltage_control.y0 = (a1 * voltage_control.y1) - (a2 * voltage_control.y2) +
                       (b1 * voltage_control.x1) - (b2 * voltage_control.x2);

  // new value for PWM
  duty = (voltage_control.y0 / PWM_GAIN) * PERIOD_VALUE;

  return duty;
}

Is this the correct way to implement a difference equation? I'm asking because the resulting controller seems unstable, and I'd like to know if it's my implementation or something else.

If you take the Z transform on both sides of the difference equation:

$$Y(z)(1-0.695z^{-1}+0.666z^{-2})=X(z)(1.001z^{-1}-0.0010155z^{-2})$$

$$G(z) = \frac{Y(z)}{X(z)}=\frac{1.001z^{-1}-0.0010155z^{-2}}{1-0.695z^{-1}+0.666z^{-2}}=\frac{2.002z^{-1}-0.002031z^{-2}}{2-1.39z^{-1}+1.332z^{-2}}$$

It does not match your transfer function. I think you can work backwards to get the right difference equation now (note the x coefficients go on the numerator and the y's on the denominator). The system you implemented has a pole outside the unit circle so it's unstable.


Finding pitch from a wideband spectrogram We can use the narrowband spectrogram for finding the pitch of an utterance. But is it possible to find the pitch using the wideband spectrogram? If so, how would that be done?

If your window length is shorter than the pitch period of a voiced utterance, the spectrogram will not be able to capture the fundamental frequency. This is the problem, yes?

But evidence of the fundamental frequency will still be available in the spectrogram because spoken voice has plenty of harmonic content. If you can identify the upper harmonics, showing up as horizontal stripes in the spectrogram, and find the size of the gap between them, this will be the fundamental pitch. Say your FFT resolution is very coarse under 200 hz, but you find a peak at 240 hz, 360 hz, and 480 hz. Well, this is a dead giveaway that the fundamental was 120 hz.

Moreover: The glottal impulse (vocal folds slamming together) is like a broadband impulse. You will get periodic vertical stripes along the spectrogram if the pitch is low enough, and your FFT window short enough.  The time difference between these stripes is the pitch period of the voice, so this is another clue you can use to determine the pitch even when the fundamental is ambiguous.


How do convolution matrices work? How do those matrices work? Do I need to multiple every single pixel? How about the upperleft, upperright, bottomleft and bottomleft pixels where there's no surrounding pixel? And does the matrix work from left to right and from up to bottom or from up to bottom first and then left to right?

Why does this kernel (Edge enhance):

turns into this image:

?

The convolution process works by multiplying every single pixel on the image with the matrix/kernel. And for the outer pixel in top, bottom, right and left there are several ways to do it:

You can ignore it, so you can start from the second line if the kernel is 3x3.
You can put a 0 value as the outer border.
You can copy the closest neighbor pixel value in as the outer border.
You can pretend the image is circular by using the opposite pixel. Top with bottom and left with right border.

As far as I know the matrix goes from left to right direction. And for that matrix (the edge enhancement), it gets pixel that have huge difference with the left neighbor. The pixel is multiplied with 1 and the left neighbor is multiplied with -1. If the two pixels have huge difference, the result will be near to white (1 or 255) but if the two pixels have small difference, the result will be near black (0).

I hope this answers your questions.


How to Do Convolution a 2D Signal (Image) and Result of Convolution of Two 1D Filters? Assume a 2-D signal (i.e., some image). Load image and assume it to be signal x. Next assume that instead of having a 2-D filter you have two one D filters

$$f_1[n] =\begin{bmatrix}0.25 && 0.5 && 0.25\end{bmatrix}$$ and
$$f_2[n]=\begin{bmatrix}0.25 \\ 0.5 \\ 0.25\end{bmatrix}$$
Assume that the convolution

$$f_1[n]* f_2[n] = f_3[n]= \begin{bmatrix}0.0625 && 0.125 && 0.0625\\ 0.125 && 0.25 && 0.125\\0.0625 && 0.125 && 0.0625\end{bmatrix}$$

Using this information and output at each stage verify that Associative property holds.

My code:

f3 = [0.0625 0.125 0.0625; 0.125 0.25 0.125; 0.0625 0.125 0.0625]
x=imread('img.jpg')
conv2('x, f3)

But this gives error saying x is not a vector. How do I fix this?

You should show that conv2(f3,x) can be implemented in a separable way. In other words, computing x1 = conv2(f1,x), then applying f2 on the results, x12 = conv2(f2,x1) gives the same result.  Or f12 = conv2(f2,f1) then x12 = conv2(f12,x).

You should be aware of the correct "dimensions" for f1 and f2: transpose of each other, so that one works on "rows" and the other on columns.

This is shown in the following Matlab code:

dataImage = zeros(25,25);
% Create a delta-like image
dataImage(13,13)=1;
% Create a random-kernel image
dataImage(11:13,11:13)=rand(3,3);
imagesc(dataImage);colormap gray

f1 = [0.25 0.5 0.25]'; f2 = [0.25 0.5 0.25];
disp(f1*f2);

f3 = [0.0625 0.125 0.0625; 0.125 0.25 0.125; 0.0625 0.125 0.0625];
disp(f3);

dataImagef1f2 = conv2(f2,conv2(f1,dataImage));

dataImagef1f2 = conv2(conv2(f2,f1),dataImage);

dataImagef3 = conv2(f3,dataImage);

subplot(2,2,1)
imagesc(dataImage);
xlabel('data')

subplot(2,2,2)
imagesc(dataImagef1f2);
xlabel('(data*f1)*f2')

subplot(2,2,3)
imagesc(dataImagef1f2);
xlabel('(f1*f2)*data')

subplot(2,2,4)
imagesc(dataImagef3);
xlabel('f3')

Suggestions for tutorials:

The Matlab Image Processing toolbox guide (900+ pages)
An Introduction to Digital Image Processing with Matlab


what is the theory behind matlab butter function? I would like to know the theory behind matlab butter function. Is there any book which gives clear idea on how it works. In future I would like to calculate the coefficients manually.

The MATLAB function that you refer to uses this knowledge to calculate the coefficients

Briefly, the function you refer to goes something like this:

The transfer functions for analog Butterworth filters of any order are already known analytically. It selects the normalized transfer function polynomial for the order you select.
The normalized polynomial is scaled appropriately according to the desired cutoff frequency, giving an $s$-domain transfer function that has the response that you want.
The $s$-domain prototype filter is then mapped to a discrete-time approximation using the bilinear transform. This is the most commonly-used method of mapping a continuous-time filter to a discrete-time filter.
The result of the bilinear transform is a transfer function in the $z$-domain, from which the filter coefficients are easily found by inspection.


Convolution of h(n) & x(n) can anyone please help me with finding Convolution for following problem .....

1) $h(n) = \delta[\sin(1 + |n|)]\quad \text{and}\quad x(n) = \sin(n^2)$

2) $h(n) = \delta[\sin(2\pi n/N)] \quad \text{and}\quad x(n) = n^2$ where $N$ is odd.

HINT:

I assume that $n$ is an integer and that we're talking about discrete-time signals. You need to simplify the expressions for $h(n)$ by noting that $\delta(n)=1$ for $n=0$ and zero otherwise. So, for which values of $n$ does $\sin(1+|n|)$ equal zero, and for which values of $n$ does $\sin(2\pi n/N)$ equal zero?


How a diode-detector affects to Frequency Response? I have a question related to crystal radio. I can't understand how a crystal detector (diode) shifts the spectrum from the carrier frequency to zero.

The antenna receives bunch of signals and the LC circuit passes the most of the frequencies to the ground excepting the resonance frequency (and the near frequencies). So it's the signal before detector:

The diode passes the positive half-waves to the capacitor/buzzer. So it's the signal after detector:

The capacitor just smooths the high-frequency ripples and the buzzer emits the original demodulated low-frequency signal.
Let's see how it works in the frequency domain.
The antenna receives signals:

The LC circuit extracts the interesting bandwidth:

The non-linear detector as i can understand should shift the bandwidth to zero (and add some other harmonics):

The capacitor/buzzer pair removes the high-frequency harmonics:

The non-linear detector thing i can't understand. How does cutting the half-wave affect to the spectrum? Why the spectrum is shifted down?
I understand how it works in the time-domain. My question is how does it work in the frequency-domain?

I tried to reproduce such behavior in the sound editor. I've generated the sine at 1000 HZ and then modulate it by the 5000 HZ carrier:

So we have 2 peaks at 4kHz (negative frequency shift) and 6kHz. The next step is the detection. So i've reproduce it by using clipping but the result spectrum isn't the same as i expect:

As you can see there is NO peak at 1kHz at all so the low-pass filter is pointless here.
Why here i've got such result? Where am i wrong?
Thanks in advance!

But what does detector in the frequency domain?

Something really complicated.  We get trained to think about signal processing almost exclusively in the frequency domain, and we forget that the frequency domain stuff was only invented as a mathematical trick to make it easier to understand what's going on in the time domain.  If it's easier to do in the time domain, do it there.
In fact, envelope detection is really horrid from a frequency domain point of view, because not only are there nonlinearities in play, but there's nonlinearities and that perfectly linear capacitor that are interacting.  So generally if you want to analyze AM radio demodulation you try to find some equivalent operation that is linear, or at least less horridly nonlinear.
In this case, just considering the demodulation to be clipping isn't bad.  Figure that the clipping operation can be denoted as $y(t) = f\left( x(t) \right )$, where $f(x)$ is (thankfully) memoryless.  So you can assert* that it has a Taylor's expansion $$f(x) = \sum_{n=0}^\infty h_n x^n$$ where the $h_n$ are the terms of the expansion**
If we feed a sum of sinusoids through $f(x)$, the above expansion tells us that we're going to get just about every possible distortion product out of the thing.  We'll get stuff at the sum of all the sinusoid frequencies, at every possible combination of sum and difference frequencies, and at ever possible combination of the sum and difference of every possible harmonic of both sinusoids.
At this point, I'm just going to claim without proof that the element at the difference between the carrier and the sideband frequencies will be prevalent -- if the carrier and the modulating signal are widely spaced apart.  AM radio works because you're transmitting at 500kHz or more, and modulating the carrier with a signal that only occupies 0Hz to 3000Hz or so.
I'm not going to try to prove this in more detail because the math gets crazy in the frequency domain -- and because you can so clearly understand the operation in the time domain.

Why the example in audio editor doesn't work?

First, because you just multiplied $\cos \omega_c t$ with $\cos \omega_s t$  For AM modulation, you want something like $\left(\cos \omega_c t\right ) \frac{s(t) + 1}{2}$ -- you multiply the carrier by something that never falls below zero.
Second, because your carrier and your modulating tone are so closely spaced (although the signal may be recognizable).
Here's some plots.

Top, a picture of one 1kHz cycle's worth of your signal, as you've defined it.  Because you didn't shift your modulating tone, you don't have enough carrier to form an envelope.

My version with your frequencies, your modulating tone shifted by 1 and reduced by half.

My version with a carrier frequency of 20kHz, and the modulating tone shifted.

* You can assert this, but you'd be lying, because it has derivatives that go to infinity -- but you can make up something arbitrarily close to a clipping function that does have a Taylor's expansion -- just assume I've done that.
** And I use $h$ in honor of the horrid things I'm doing to the math...


Band limited signal: Infinite time duration, time limited signal : infinite bandwidth In signal processing, there is a time-frequency duality. Band-limited signals have infinite time duration and time-limited signals have infinite bandwidth.

Where can I find a proof of this and is there some intuitive explanation for this?

I am not sure if it has to do something with fouriers uncertainty principle?

Here is an argument that doesn't require the uncertainty principle.

Assume $y(t)$ is a signal that is not identical to zero, with Fourier transform $Y(f)$. Let $s(t) = \text{rect}(t/T)$ a rectangular pulse of duration $T$ and Fourier transform $S(f) = |T| \text{sinc}(T f)$.

Then, $x(t) = y(t)s(t)$ is a time-limited signal with Fourier transform $X(f) = Y(f) \ast S(f)$. Since $S(f)$ has infinite support, so does $X(f)$.

Since any time-limited signal $x(t)$ can be defined in this way for a suitable choice of $y(t)$ and $s(t)$ (and a possible time-shift that doesn't affect the result), then we can say that all such signals have infinite bandwidth.

The same argument can be applied for the inverse Fourier transform: band-limiting a signal involves multiplication by a rectangular pulse in the frequency domain, or a convolution with a sinc in the time domain, resulting in an infinite-duration signal.


Simplest way to implement FIR filter in C/C++ What's the simplest way to implement a FIR filter in C/C++ ?

(I don't want to have to use a third-party software like Matlab or Octave in order to compute the coefficiens because I have to use different filters on thousand of files, with different frequencies for each file.)

Thus, is there a ready-to-use C/C++ DSP library with functions like file.FIR_NotchFilter(frequency,Qfactor) ? or something similar ?
(commercial or opensource)

Thanks a lot in advance.

PS1 : I see that there is http://aquila-dsp.org, but there are no binary releases yet, and for older version, it seems difficult to use it

PS2 : A really ready-to-use solution exists for IIR filters : https://github.com/vinniefalco/DSPFilters

// Create a Chebyshev type I Band Stop filter of order 3
Dsp::SimpleFilter <Dsp::ChebyshevI::BandStop <3>, 2> f;
f.setup (3,    // order
     44100,// sample rate
     4000, // center frequency
     880,  // band width
     1);   // ripple dB
f.process (numSamples, arrayOfChannels);

Have you considered using the Parks McClellan algorithm to generate your FIR filters. The source code for it is available on several sites in Fortran or C.

The original Fortran code is available on Wikipedia.

Here are two sources for C code

Github/scipy
Iowa Hills

You said you were trying to understand filters, so let me explain that the Parks McClellan algorithm is a classic. It was developed in the early 70's, and is still considered to be one of the best FIR design algorithms available.


Why is Error Vector Magnitude (EVM) measured as an RMS value? I do understand Voltage/Current measurement of a Sine Wave is in RMS, as average yields a zero, and RMS gives a true 'DC equivalent' value. But I do not understand why an EVM measurement could be in RMS, and not simply Average.

Is it because, similar to a sine wave, as symbols are mirrored in the IQ constellation plane, their individual EVM values when averaged yield a zero as well?

You've got it – zero mean noise still has power. That power happens to be its variance. RMS of such signals thus happens to be the noise amplitude's standard deviation – and give a sensible number to assess link quality and error probabilities.


Simple way to detect continuous lines I have some images of a background, and of steam that I am subtracting to extract just the steam. Unfortunately, because the camera conditions are not perfect the differenced image is always a few pixels off generating sharp lines where features in the background don't match up pixel-by-pixel like so

I'm interested in segmenting out the steam (foreground noise), from the background. I think if I can automatically remove the continuous lines from this image it would be sufficient for my purposes.

Unfortunately it's been almost ten years since my image analysis work and I don't remember what you do to detect continuous lines. What are some simple well preforming algorithms to do this?

I have two methods to suggest:
1. Use the erosion (morphological operator) (you can also use other operators like closing/opening). But i think it will destroy the steam data also...
2. A better solution is to detect the lines using hough transform. Google it and you'll find tons of info about it (if your'e using MATLAB, there is a guide how to do it in their website)

Have fun...


Autocorrelation - cyclic vs acyclic In terms of signal analysis (especially frequency decomposition) would it be bad to use cyclic autocorrelation (AC) instead of acyclic/linear AC?

My use case:
I'm trying to characterize events in a sensor stream by periodicity/frequency. At the moment, I'm using DFT spectrum (spectral density) and AC of sliding windows.
On the one handside it's algorithmically way cheaper to calculate the cyclic AC via FFT than acyclic (pad with zeroes and do double length FFT).
However, the cyclic AC looks similar (though not equal) to the acyclic for the first N/2 lags and is self-symmetric - thus at most half the information of acyclic AC can be contained.

My motivation:
If the cyclic AC would be sufficient for detection, that would mean better performance compared to linear AC (in terms of memory usage and computation time). I then could also calculate AC inplace (better parallelization potential).

Rephrased question:
Is there any rule when cyclic AC should NOT be used?
Or perhaps, are there known cases, when cyclic AC might actually be better than linear AC?

Thanks in advance :-)

If the cyclic autocorrelation output has the properties you need in your application (e.g., it's peak-to-average ratio is high enough for picking out the peaks), then it's probably okay to substitute the cyclical correlation for the linear one. For example, consider a length 7 maximal-length sequence:
$$
\{+1, -1, -1, +1, +1, +1, -1\}
$$
The cyclical and acyclical autocorrelation functions are plotted below.

Notice that both functions have a peak value of 7. The acyclical autocorrelation function has a largest off-peak magnitude of 3, whereas the cyclical autocorrelation function has a largest off-peak magnitude of 1. The peak-to-off-peak amplitude ratios are $7$ and $7/3$ for cyclical vs. acyclical, respectively. This means that the cyclical autocorrelation function is basically three times more peaky, making it much easier to detect. If the extra peak-to-off-peak amplitude ratio is not necessary for your application, it would be okay to correlate using linear correlation instead.

If instead, you have something whose acyclical autocorrelation properties are optimal (such as a Barker code), and you would like to calculate that acyclical correlation but don't want to use the computationally inefficient direct form of convolution, there are algorithms that calculate the acyclical autocorrelation using the computational savings of FFTs. Two of these I know of are called the overlap-add and overlap-save methods. They work on blocks of the input data at a time, and using FFTs of the blocks of data and multiplication in the transformed domain, calculate the acyclical correlations/convolution of the data against your given sequence.


How to estimate the modulation transfer function of images? I'd like some way to estimate the modulation transfer function in images.  This is the observed drop in amplitude vs increase in frequency of spatial details - nice illustration here:

(from https://www.imatest.com/docs/sharpness/)
In my case, I have a large collection of digital images of the same kind of thing; not literally the same object, but a fairly uniform set of objects with consistent overall properties (eg amount of spatial detail).  The imaging system is also exactly the same, but (because this is a fairly challening setup) the image quality varies considerably from image to image.  One of the ways in which it varies is that some parts of some images are seriously out of focus.  These are large images, around 10MP, and there is a lot of them (tens of thousands) so the statistical properties over the whole set can be determined with some confidence.  I think the typical image can resolve spatial detail on the scale of a couple of pixels (let's say 2-3 times the pixel pitch) but in some examples it looks like what it would look like if that area had a gaussian blur with radius ~10 pixels.  The images contain some sharp edges (features which are expected to be a fairly sharp edge in reality - they don't always look very sharp in the images) but no particular area is guaranteed to have an edge.
I'd like to be able to calculate an approximation to the real modulation transfer function just from the real images (or parts of images), without having something like a test chart in the image.  I think the consistency of the image statistics should make that possible in principle, but I'm not sure how to approach it.
How can this be done?

The "classic" way of estimating blurriness is (as far as I know) the difference of Gaussians approach: you filter your image with Gaussian blur filters of different variance. Then, you at which level of blurriness the most info is lost (i.e. where the largest difference between the last level of blurriness and the next level is).
If your image content is reliably high-frequency to begin with, a simple high pass filter followed by energy detection on regions (i.e. parts) of the image can be used to.
In the end, if you combine the idea of mapping linear operations (Gaussian blur, high pass filters) to spatially constrained regions, you end up with wavelets. This might be a very promising approach if you need to analyze properties of images later on (think, as a simple example, something like "images which are sharp in about ¼ of the area", or "images that contain both very blurred and very in-focus regions" or so).


Why do the two methods give different answers for the Fourier transform of $Y = \cos(\omega_0 t + \phi)$? Why do the following two methods give different answers (or are they the same) for the Fourier transform of $Y = \cos(\omega_0 t + \phi)$, with respect to $t \to \omega$ ?

They are the same.

Using the Dirac impulse sifting property $f(x) \delta(x-a) = f(a)\delta(x-a)$ you can also verify that the second method produces the same output as the first method.


RMS in frequency domain using Plancherel theorem I have an acceleration measurement in x, y and z direction, i.e. three vectors each of length 3000.
$x = [x_{1}, x_{2}, x_{3},...,x_{3000}]$
$y = [y_{1}, y_{2}, y_{3},...,y_{3000}]$
$z = [z_{1}, z_{2}, z_{3},...,z_{3000}]$
Currently I'm taking the RMS of x, y and z then I'm calculating the FFT. This produces a vector of length 3000.
$RMS = [\sqrt{x_{1}^2 + y_{1}^2 + z_{1}^2}, \sqrt{x_{2}^2 + y_{2}^2 + z_{2}^2},...,\sqrt{x_{3000}^2 + y_{3000}^2 + z_{3000}^2}]$
I've now read that it is also possible to apply the RMS in frequency domain then using the Plancherel theorem which states that:
$\sum_{n=0}^{N-1}|x_{n}|^2 = \frac{1}{N}\sum_{k=0}^{N-1}|X_{k}|^2$
if $X_{k}$ is the FFT of $x_{n}$.
As far as I understand it this would produce a single output for my three vectors of length 3000. Is that correct? Is there a possiblity to compute the RMS in frequency domain and getting a vector of length 3000?

The array $[\sqrt{x_{1}^2 + y_{1}^2 + z_{1}^2}, \sqrt{x_{2}^2 + y_{2}^2 + z_{2}^2},...,\sqrt{x_{3000}^2 + y_{3000}^2 + z_{3000}^2}]$ is a vector of the Euclidean lengths of 3000 3D-space vectors.
The designation RMS (root mean square) is usually used for the square root of the arithmetic mean of the squares of serial values of the same kind. It is not used to designate the value $\sqrt{x^2 + y^2 + z^2}$, which is called the space vector length. The values used for RMS calculation can be received from a series of measurements of the same observable. The components of a 3D vector are not three instances of measurement of an observable; in fact, these vector components are viewed as three observables: the x component, the y component, and the z component.
You have 3000 samples from a series of measurements of acceleration, each instance of measurement produces three vector components of a 3D acceleration vector. When calculating RMS, you iterate over individual items which you select from a vector of measurement samples -- selecting one sample of a total of 3000 samples, squaring the selected item value, summing up these squares in the course of iteration -- then you divide the total sum by the number of samples (3000 here), and finally calculate the square root. You calculate three RMS values: an RMS for the x component of the selected acceleration vector, an RMS for the y component of the same vector, and an RMS for the z component of the same selected vector, like this:
$$
RMS_x = \sqrt{{1\over N}\sum_{n=0}^{N-1}x_{n}^2}; \,
RMS_y = \sqrt{{1\over N}\sum_{n=0}^{N-1}y_{n}^2}; \,
RMS_z = \sqrt{{1\over N}\sum_{n=0}^{N-1}z_{n}^2}
$$
These $RMS_{x,y,z}$ components are components of a root-mean-squared acceleration vector, $\overrightarrow {RMS} = [RMS_x, RMS_y, RMS_z]$
FFT and the Plancherel theorem may or may not be required for processing of your data, dependent on what processing you are going to do on the measurement data.
Notice the usage of the word vector for two different entities: the vector as an element of generalized vector space (as a set of subsequent samples from a waveform) versus the Euclidean vector (a directed segment). See https://en.wikipedia.org/wiki/Vector_(mathematics_and_physics).
It may happen that you need to compute the RMS of acceleration vector lengths, irrespective to acceleration direction in space. If this is the case, the formula is
$$
RMS_{scalar} = \sqrt{{1\over N}\sum_{n=0}^{N-1}(x_{n}^2 + y_{n}^2 + z_{n}^2)}
$$
subscript scalar means that, unlike $\overrightarrow {RMS}$, $RMS_{scalar}$ it is not a vector. Notice also that $RMS_{scalar} = \sqrt{RMS_x^2 + RMS_y^2 + RMS_z^2}$


Apply FIR filters to Extract Stego from WAV data I appreciate anyone that takes a moment to help me with this problem. I've been banging my head against the keyboard for a while, searching forums and DSP tutorials and I can't figure this problem out.

Here's a couple links I found very helpful from this forum on the subject:
http://www.dspguide.com/
http://wiki.scipy.org/Cookbook/ApplyFIRFilter

So I have a WAV file with a 28 sec song/sound that I am trying to convert to frequency range, apply an FIR filter to isolate the 14.5 - 16kHz range. The WAV has a sampling rate of 32 kHz and contains an intermittent buzz that becomes clearer when I drop out all other frequencies besides the 14.5 -16kHz bands.

I have written this sample of python code, and believe I am close to isolating the signal. Its hard to tell.

# read data from file name fn, using low/high for bandpass
# low and high frequency cutoffs, and a
def filter_frequencies(fn,low=0,high=0):
    data = wave.open(fn,'rb'); rate = data.getframerate()
    ntap = (rate / 4) + 1; nyq = rate / 2
    low,high = low/nyq,high/nyq
    chunks = int(data.getnframes()/(rate*4))
    ns = b''
    if low and high == 0: # lowpass
        bp = low; pz = True
    elif high and low == 0: # highpass
        bp = high; pz = False
    else: # bandpass
        bp = [low,high]; pz = False

    b = firwin(ntap, bp, pass_zero=pz, window='blackman', nyq=nyq)
    y = b''
    for num in range(chunks):
        # read 4 secs of frames from wav file
        seg = np.fromstring(data.readframes(rate*4),dtype=np.int8)
        #seg = np.fromstring(data.readframes(data.getnframes()),dtype=np.int8)
        # convert frame amplitudes to their frequency values
        freq = np.fft.rfft(seg)
        # perform lowpass, highpass or bandpass filter
        cr = fftconvolve(freq,b,mode='valid')
        # converts frequencies to binary string
        y += cr.tostring()
    data.close()
    return y

I am a noob to DSP, so again I would really appreciate your help on if I am doing the FFT to inverse FFT correctly.

Below is the code I am using to transform the output from the above function back into WAV bytes format.

def frequency_to_bytes(freq):
    # invert the Fast Fourier Transform
    nf = np.fft.irfft(freq)
    # convert frequencies back to integer values
    bs = np.ravel(nf).astype(np.int8)
    return bs.tostring()

P.s. if you have seen this problem before, or know what challenge I am trying to solve, please no spoilers. I have been working on this for a week solid.

edit: added some background links I have been using as reference.

I never worked with Python or SciPy, but I don't see the sense in the actions you apparently perform in your code.
firwin provides you with the time-domain filter coefficients, which are to be convolved with the time-domain signal you wish to filter.
What you do in your code, is convolve the time-domain filter coefficients with the frequency-domain signal - which makes absolutely no sense.

What is the point of converting the signal to freq-domain (doing FFT)? Is the stego analysis/decoding performed in time- or frequency-domain?

If your processing chain requires conversion of the signal to freq-domain, for some purpose, then you don't really need the FIR filter, since you can get rid of the undesired frequency elements once you have the signal in the frequency domain, by simply zeroing-out the undesired frequency bins of the signal.

If the only reason for calculating FFT of the signal is to filter it - then there is no need to do it, as you can simply convolve the time-domain filter coefficients with the time-domain signal, as I explained.


How to plot BER Vs SNR of a received signal I am using USRP to transmit and received MIMO frames similar to Wi-Fi standard that I created on my own. After receiving the signal, it is processed as described in the standard.

While I plot the QPSK constellation points after frequency and phase correction, I have a plot (shown above) which is noisy. Now I want to plot (BER Vs SNR).  I know that we can calculate BER as: (Error Bits/Total No. of bits). But I am not sure how to calculate SNR from these complex symbols.

What I understand from online materials is that SNR is the noise variance that is estimated from around the constellation points (maybe some distance/deviation measurement form the true point. Or its is given by Eb/N0, where Eb is the energy of bits received and N0 is the noise. The question might be very silly, what I do not understand is how I should find Eb and No from these complex received values that is plotted on a constellation diagram. I am doing all the coding in python.

From the plot it seems the transmitted QAM points are $(0.5, 0.5),(-0.5,0.5),(-0.5,-0.5),(0.5,-0.5)$, and if you can model your comm system as just AWGN channel, then your signal amplitude will just be $\sqrt{0.25+0.25} = \sqrt{0.5}$ and hence signal power will be $0.5$. In doing so, I have assumed all 4 QAM symbols equiprobable.

You can get your noise power by averaging all squared Euclidean distances of received points from their true QAM points.

From this you can get received SNR.


When is the Fourier transform of a signal periodic? I mean not the time-domain signal being periodic, but the Fourier transform being periodic.

I'm new to this exchange and I'm not sure how mathy you all get. I think the answer below is cool because it shows that in some sense the continuous-time Fourier transform is never periodic but that in another sense there are lots of ways to get periodic transforms.

For the continuous-time Fourier transform on $\mathbb{R}$, both CMDoolittle's and Robert Bristow-Johnson's answers are correct, though for slightly-different definitions of "Fourier Transform". In the usual way we talk about the continuous Fourier transform, transforms are never periodic. If we're willing to transform distributions (of which the Dirac delta is an example), we can have transforms which are periodic.

Usually, the theory of the continuous Fourier transform is developed on $L^2(\mathbb{R})$. The Fourier transform is a bijection of $L^2(\mathbb{R})$ back onto itself; this means that $L^2(\mathbb{R})$ is also the space of all possible Fourier transforms. However, the zero function is the only periodic function in $L^2(\mathbb{R})$, so we can conclude that continuous Fourier transforms of non-zero functions are never periodic. This is CMDoolittle's answer.

However, we also have a Fourier theory for distributions. Specifically, we define distributions to be the set of linear functionals on the space of so-called test functions. Test functions, in turn, are $C^\infty$ functions which are required to decay rapidly (alternately, sometimes test functions are $C^\infty$ functions with compact support). Some of these distributions can be identified with periodic functions on $\mathbb{R}$. For example, if $f$ is a test function then the mapping
$$ S : f \mapsto \int_{-\infty}^\infty \sin(t) f(t)\ dt$$
is a distribution which can be identified with the sine function. As we can define the Fourier transform of a distribution, so we have a way of defining the Fourier transform of periodic functions. As it turns out, the Fourier transform is also a bijection of the space of distributions back onto itself. Going back to the distribution $S$, this means that we can write $S$ as the Fourier transform of some other distribution $T$: $S = \hat{T}$; in other words, we conclude that there has to be a distribution $T$ whose transform is the distribution $S$ which is identified with the periodic function sine. Clearly, we could pull this same trick with, say, any bounded periodic function. Ultimately then, we can extend Robert Bristow-Johnson's answer: if we're talking about the Fourier transform of distributions, there are lots of distributions whose Fourier tranforms are periodic.

For a discussion of the Fourier theory of distributions, see Chapter 7 of Rudin's Functional Analysis.


Why is my time domain interpolation via zero-padding in frequency domain wrong? Since the process can be applied in either domain to increase the sampling rate in the other domain, I am trying to apply zero-padding in frequency space to recover a 'cleaner' interpolated signal in temporal space.
To do so, I insert zero-valued frequencies in the spectrum at the location of higher frequencies, which is a common practice.
However I don't seem to recover the original signal very well (in black below) after zero-padding (in red).
import numpy as np
import matplotlib.pyplot as plt

# odd dimension for simplicity
n    = 19
npad = 99

x    = np.linspace(0.,4.*np.pi,n)
xpad = np.linspace(0.,4.*np.pi,npad)

f = np.cos(x) + 1j*np.sin(x)

f_fwd = np.fft.fft(f)

f_fwd_pad = np.zeros(npad,dtype=complex)

h = (n-1)//2
f_fwd_pad[0:h+1]   = f_fwd[0:h+1]
f_fwd_pad[npad-h:] = f_fwd[h+1:]

f_interpolated = np.fft.ifft(f_fwd_pad)*npad/n

fig, ax = plt.subplots(1,2)

ax[0].plot(x,np.real(f),linestyle=None,marker='x',color='k')
ax[0].plot(xpad,np.real(f_interpolated),color='r')

ax[1].plot(x,np.imag(f),linestyle=None,marker='x',color='k')
ax[1].plot(xpad,np.imag(f_interpolated),color='r')

Is that result expected? Is there some fundamental understanding that I am missing?

The answer above is correct. Just to clarify a bit further, using x = np.linspace(0,10,5) will produce 5 numbers from 0 to 10 inclusively
np.linspace(0,10,5)
array([ 0. ,  2.5,  5. ,  7.5, 10. ])

You don't want the last number because in your example the last number is the first number of the next period. A correct implementation would be:
periods = 4.*np.pi
x    = np.arange(0., periods, periods/n)
xpad = np.arange(0., periods, periods/npad)


ICY feature detector and edge detector in opencv I need to trace the internal feature boundaries of blobs (blood cells) in OpenCV to train a classifier.

As a good reference from ICY in java, the following sequence of operations gives good, workable edge contours.

Edge Detector( Feature Detector( CLAHE (Original Image) ) )

[UPDATE1 @Amitay] In the above image, the first reddish blob is the original image, the second one is the result of the Feature Detector operation from ICY and the third one is the result of applying the Edge Detector on the detector features.

Attempting this in OpenCV using the LabEdgeDetector.cpp ( by Austin Walters) followed by Canny Edge detector, I get the following - not so smooth results. There are lots of X and Y directional artifacts throwing the edge detection off.

From what I can see in the java source, the ICY feature detector plugin uses steerable filters for its operation (steerabledetector.java).

Are there equivalent featuredetector algos in OpenCV?  If so please list them.
How can I improve on my results? Please guide me on any other sequence of existing OpenCV operations that can produce similar results to those from ICY above.

p.s I dont have enough reputaions to post more links and references.

I would suggest to smooth the image before you apply the edge detection.
Apply bilateral filter on the image before the edge detection( cv2.biltralfilter).


Sin wave spectrum shows multiple frequency I am very new to DSP, so excuse my ignorance.

Using Audacity, I generated a sine wave of:

frequency 440Hz
amplitude: 0.8
duration: 1 minute

Now when i plot the spectrum for this sine wave, I get an image like this:

This is really surprising for me.

I was expecting a single surge at 440 Hz for the pure tone. But this is entirely different.
Sure this is not anything to do with noise. Right ?

My question

How can the spectrum of a pure tone have so many frequency components in the spectrum.

What obvious am I missing ?

An FFT only results in a single component (plus its complex conjugate and some tiny numerical noise) if the single pure real sinusoid input is of a frequency whose period is an exact integer submultiple of the FFT size.  A rectangular window makes any "leakage" into other non-adjacent bins even worse.

But the length of the FFT you have specified is 512, which is not an exact multiple of the period of your signal (which is 44100/440 = 100.227273).


If video uses color bars to test and align color values among TV sets, what could be the equivalent in Audio? Are there any sample Audio files, encoded in the simplest format known, that can mathematically explain or illustrate the fundamentals of (audio) signal processing that might also relate to known equations of physics pertaining to audio?
I am trying to teach a robot AI system how to sing.

If you want the AI to relate to how we perceive sound, maybe something like

a linear chirp from 20Hz-20kHz
white noise in the same bandwidth
a single sine wave at 1kHz (roughly our greatest sensitivity)
a harmonic series from 55hz doubling up to 14.08kHz (all A notes to illustrate octaves)
other harmonic series and/or progression to illustrate important musical intervals (major and minor chords, pentatonic scale)
sound level variation of one or more of the above between 0dB SPL and some appropriate limit to establish dynamics

In terms of storage, a set of vectors is probably the best place to start. You can make one long vector with your test pattern and load it into and out of a WAV file I suppose.


Analysing mp3 quality with spectrum, how to interpret the results? When comparing spectrum images between two different audio sources in order to determine which has better quality, I sometimes come to an impasse. For example, consider the following images from two 320kbps mp3 files:

Optical difference of RGB values:

(here is the difference between waveform-synced versions of the above examples, but it's quite the same)

The second version has a cutoff at 20 KHz, while the first contains some audio information above that frequency, albeit very low in intensity, and probably inaudible.

On the other hand, this extra audio portion at the top of the plot appears to be noise, since it's composed by an almost continuous "blurry" zone. In addition, looking carefully at the two images, some subtle differences can be observed at the lower frequencies.

I'm inclined to think that the second version is better quality, theoretically, since the lower cutoff limit may allow for the encoder to make better use of the available bandwidth, so that it can more acuratelly represent the lower frequencies. That aside, I really don't know how to interpreted those subtle differences between the two images in order to make an objective decision.

So, when comparing spectrum images, what specific details can help determine which encoder did a better job?

At those relatively high rates, the possibility (for a human) to hear a difference depends a lot on the device it is played on. My first suggestion would be to rely on a panel of testers (possibly bats, for the high frequencies), use the best audio system you can, and compute a Mean Opinion Score.

I do agree that the frequency limitation around 20 kHz might seem in disadvantage of the second image, but as you said, low-freq bits can then be better allocated (possibly some).

Similarly as @Marcus Müller, I suspect the high-freq band in the first to be spurious. Especially since it really looks symmetric in the time-frequency plane (modulo the contamination with the real audio), and the aliasing hypothesis sounds nice here.

As said in a comment, you can subtract one decompressed sound from the other. Hopefully, by a careful listening, you could check whether the difference is meaningful, or not.

So, I'll, too,  go for the second for the best.