Optimal amplitude of an $m$-bit sinusoid A continuous-time sinusoid of zero-to-peak real amplitude $A \le 2^{m-1}-0.5$ (e.g., for $m=16$, $A \le 32767.5$) is quantized to $m$-bit resolution by rounding it to the nearest integer (Fig. 1). What is the optimal amplitude for different $m \le 24$?

Figure 1. For $m = 4$, a sinusoid of the maximum allowed amplitude $A = 7.5$ (blue solid line), its quantization (orange dashed line), and the most negative and positive 4-bit two's complement signed integers -8 and 7 (yellow).

An answer is desired for each of two alternative definitions of optimality:

maximize signal-to-noise ratio (SNR), where the signal is the original sinusoid and noise is the quantization error, and
minimize total harmonic distortion (THD), which is normalized against the amplitude of the fundamental frequency in the Fourier series of the quantized sinusoid.

The phase of the sinusoid does not matter: A phase shift of a sinusoid is equivalent to shifting it in time, which results in a time shift of both the quantized sinusoid and the quantization error. The power spectrum is invariant to time shifts. We choose to work with sinusoid $A\cos(x)$.

Optimality by def. 1

Equivalently to maximizing signal-to-noise ratio (SNR), we can minimize the root mean square quantization error relative to the root mean square $\sqrt{1/2}A$ of the sinusoid, $\mathrm{RRMSE}$ $=$ $\sqrt{1/\mathrm{SNR}}$. It is enough to do the analysis in the first quarter-period of the cosine wave, because the rest of both the cosine wave and the quantization error are identical to their first quarter-period up to a sign flip and/or time reversal. Let $x$ be in the first quarter of the cosine wave, $0 < x < \pi/2$. Along the lines of Eq. 2 of my answer to a related question, the quantized cosine wave attains an integer value $k \in 0\ldots\mathrm{round}(A)$ when:

$$\begin{align}0 < x < \mathrm{acos}\left(\frac{\mathrm{round}(A) - 0.5}{A}\right),\quad &\text{if }k = \mathrm{round}(A),\\
\mathrm{acos}\left(\frac{k + 0.5}{A}\right) < x < \mathrm{acos}\left(\frac{k - 0.5}{A}\right),\quad &\text{if }1 \le k \le \mathrm{round}(A) - 1,\\
\mathrm{acos}\left(\frac{0.5}{A}\right) < x <\frac{\pi}{2},\quad &\text{if }k = 0.\end{align}\tag{1}$$

Each value of $k$ gives to the relative mean square error $\mathrm{RelMSE}$ an additive contribution of:

$$\frac{2\,\mathrm{MSE}_k}{A^2} = \frac{2}{A^2}\frac{1}{\pi/2}\int_{x_0}^{x_1}\big(A\cos(x)-k\big)^2dx = \frac{2\sin(x_1)\cos(x_1) - 2\sin(x_0)\cos(x_0)}{\pi} + \frac{8k\big(\sin(x_0) - \sin(x_1)\big)}{\pi A} + \frac{4k^2(x_1 - x_0)}{\pi A^2} + \frac{2x_1 - 2x_0}{\pi},\tag{2}$$

where $x_0$ and $x_1$ denote, as defined separately for each $k$, the bounds $x_0 < x < x_1$ given by Eq. 1. The total relative mean square quantization error is then:

$$\begin{gather}\mathrm{RelMSE} = \frac{2\,\mathrm{MSE}}{A^2} = \sum_{k=0}^{\mathrm{round}(A)}\frac{2\,\mathrm{MSE}_k}{A^2}\\
= \frac{2\,\mathrm{asin}\left(\frac{1}{2A}\right)}{\pi} - \frac{\sqrt{4A^2 - 1}}{2\pi A^2} - \frac{(2A^2 + 4\,\mathrm{round}(A)^2)\,\mathrm{asin}\left(\frac{2\,\mathrm{round}(A) - 1}{2A}\right)}{\pi A^2} - \frac{(6\,\mathrm{round}(A) + 1)\sqrt{4A^2 - (2\,\mathrm{round}(A) - 1)^2} - 2\pi(A^2 + 2\,\mathrm{round}(A)^2)}{2\pi A^2}\\{+\,\frac{1}{\pi A^2}\sum_{k=1}^{\mathrm{round}(A)-1}}\left(\left(2A^2 + 4k^2\right)\bigg(\mathrm{asin}\left(\tfrac{2k + 1}{2A}\right) - \mathrm{asin}\left(\tfrac{2k - 1}{2A}\right)\bigg) + \frac{(6k - 1)\sqrt{4A^2 - (2k + 1)^2} - (6k + 1)\sqrt{4A^2 - (2k - 1)^2}}{2}\right),\quad A > 0.5.\end{gather}\tag{3}$$

For optimality by definition 1, the task is to find the value of $A$ that minimizes the relative root mean square quantization error RRMSE $= \sqrt{\mathrm{RelMSE}}$. under the constraint $A \le 2^{m-1}-0.5$. Eq. 3 can be evaluated in Python using the mpmath arbitrary precision library:

import mpmath as mp
def RelMSE(A):  # valid for A >= 0.5
    A = mp.mpf(A)
    return 2*mp.asin(1/(2*A))/mp.pi - mp.sqrt(4*A**2-1)/(2*mp.pi*A**2) - (2*A**2 + 4*mp.floor(A + 0.5)**2)*mp.asin((2*mp.floor(A + 0.5) - 1)/(2*A))/(mp.pi*A**2) - ((6*mp.floor(A + 0.5) + 1)*mp.sqrt(4*A**2 - (2*mp.floor(A + 0.5) - 1)**2) - 2*mp.pi*(A**2 + 2*mp.floor(A + 0.5)**2))/(2*mp.pi*A**2) + mp.nsum(lambda k: (2*A**2 + 4*k**2)*(mp.asin((2*k+1)/(2*A)) - mp.asin((2*k-1)/(2*A))) + ((6*k-1)*mp.sqrt(4*A**2 - (2*k + 1)**2) - (6*k + 1)*mp.sqrt(4*A**2 - (2*k - 1)**2))/2, [1, mp.floor(A + 0.5)-1])/(mp.pi*A**2)

RRMSE appears to have a local minimum between each pair of successive integer $A$ (Fig. 1).

Figure 1. RRMSE (blue solid line and blue squares) and its approximation $\sqrt{1/6}/A$ (orange dashed line) based on the variance $1/12$ of a uniform distribution of width 1, for various ranges of $A$.

A selection of optimal $A$ are given in Table 1 together with the resulting RRMSE, also for some other common values of $A$. At larger $m$, RRMSE reduction by the optimal choice will be marginal. The table entries, which show only digits that are equal between two computations using different precision settings, can be generated by the following Python script (continued), the execution of which took days:

def approx_optimal_A(m):
    m = mp.mpf(m)
    return 2**(m-1) - 1 + mp.mpf("0.156936321399") + mp.exp(-1.2749819017 - 0.3464088917*m)  # Eq. 8
    # return 2**(m-1) - 1  # This less informed guess gives identical results but slower convergence

def to_max_digits(f, prec_1, prec_2, max_digits):  # return the at most max_digits digits of function f that are agreed about by both precision settings
    prec = mp.mp.prec
    mp.mp.prec = prec_1
    y_prec_1 = f()
    mp.mp.prec = prec_2
    y_prec_2 = f()
    digits = max_digits
    while mp.nstr(y_prec_1, digits, strip_zeros=False) != mp.nstr(y_prec_2, digits, strip_zeros=False):  # Beware: a possible infinite loop
        digits -= 1
    return mp.nstr(y_prec_2, digits, strip_zeros=False)
    prec = mp.mp.prec

double_digits = 15  # Print at most this many digits
dB_digits = 9

for m in range(2, 25):
    optimal_A = to_max_digits(lambda: mp.findroot(lambda A: mp.diff(RelMSE, A), approx_optimal_A(m)), 80, 100, double_digits)
    RelMSE_optimal = to_max_digits(lambda: 10*mp.log10(RelMSE(mp.mpf(optimal_A))), 80, 100, dB_digits)
    RelMSE_1 = to_max_digits(lambda: 10*mp.log10(RelMSE(2**(m-1)-1)), 80, 100, dB_digits)
    RelMSE_2 = to_max_digits(lambda: 10*mp.log10(RelMSE(2**(m-1)-0.5)), 80, 100, dB_digits)
    print(str(m)+"&"+optimal_A+"&"+RelMSE_optimal+"&"+RelMSE_1+"&"+RelMSE_2+"\\\\")

Table 1. Optimal $A$ by definition 1 for different $m\le24$ and the resulting RRMSE, with the RRMSE for some common choices of $A$ listed for comparison. $m=1$ has been left out because it cannot be handled by Eq. 3 and because a single bit does not represent any positive numbers as a two's complement representation. Note that RRMSE in dB is converted to SNR in dB by a sign flip, because $\mathrm{SNR}$ $=$ $1/\mathrm{RelMSE}$ $=$ $1/\mathrm{RRMSE}^2$.

$$\begin{array}{l|l|lll}m&\text{optimal }A&\rlap{\text{RRMSE (dB)}}&&\\
&&A=\text{optimal}&2^{m-1}-1&2^{m-1}-0.5\\
\hline
2&1.26827949461530&-11.1128053&-8.92729805&-10.1764645\\
3&3.23800942121037&-18.8206937&-17.9588863&-17.8882004\\
4&7.21658597929407&-25.5167673&-25.0903048&-24.7375654\\
5&15.2007181331955&-31.8127094&-31.5775537&-31.2013629\\
6&31.1888756714257&-37.9364880&-37.7977883&-37.4726954\\
7&63.1800835394190&-43.9858910&-43.9001114&-43.6414308\\
8&127.173613625523&-50.0049518&-49.9500053&-49.7527817\\
9&255.168894736361&-56.0134743&-55.9773382&-55.8307200\\
10&511.165479188890&-62.0200022&-61.9957638&-61.8884975\\
11&1023.16302205377&-68.0278661&-68.0113689&-67.9337180\\
12&2047.16126264484&-74.0380590&-74.0267095&-73.9708981\\
13&4095.16000722516&-80.0505958&-80.0427267&-80.0028089\\
14&8191.15911371601&-86.0651409&-86.0596541&-86.0312009\\
15&16383.1584789666&-92.0812822&-92.0774410&-92.0572079\\
16&32767.1580286428&-98.0986407&-98.0959437&-98.0815798\\
17&65535.1577094659&-104.116904&-104.115007&-104.104821\\
18&131071.157483397&-110.135830&-110.134493&-110.127276\\
19&262143.157323352&-116.155234&-116.154291&-116.149181\\
20&524287.157210089&-122.174984&-122.174318&-122.170701\\
21&1048575.15712995&-128.194979&-128.194509&-128.191950\\
22&2097151.15707326&-134.215150&-134.214818&-134.213008\\
23&4194303.15703317&-140.235446&-140.23521&-140.233931\\
24&8388607.15700481&-146.255831&-146.2557&-146.25476\end{array}$$

Optimality by def. 2

A periodic function such as a quantized sinusoid has a Fourier series; it is a sum of harmonic sinusoids, that is, sinusoids of harmonic frequencies of a fundamental frequency. Harmonic sinusoids are orthogonal. Therefore, the mean square of the periodic function equals the sum of the mean squares of the harmonic sinusoids. The mean square of the sum of the non-fundamental harmonics can then be calculated by subtracting the mean square of the fundamental from the mean square of the periodic function. For a quantized cosine wave $\mathrm{round}\big(A\cos(x)\big)$ this allows to calculate total harmonic distortion (THD) as:

$$\mathrm{THD} = \sqrt{\frac{\displaystyle\mathrm{MS} - a_1^2/2}{a_1^2/2}} = \sqrt{\frac{\mathrm{MS}}{a_1^2/2} - 1},\tag{4}$$

where MS is the mean square of the quantized cosine wave, $a_1^2/2$ is the mean square of the fundamental, and $a_1$ is the coefficient of the fundamental frequency cosine in the Fourier series of $\mathrm{round}\big(A\cos(x)\big)$, calculated using Eq. 3 of my answer to a related question and simplifying to:

$$a_1 = \frac{2}{\pi A}\left(\mathrm{round}(A)\sqrt{4A^2 - \big(2\,\mathrm{round}(A) - 1\big)^2} + \sum_{k=1}^{\mathrm{round}(A)-1}k\left(\sqrt{4A^2 - (2k - 1)^2} - \sqrt{4A^2 - (2k + 1)^2}\right)\right).\tag{5}$$

Mean square of the quantized cosine wave is calculated in its first quarter-period by:

$$\mathrm{MS} = \frac{1}{\pi/2}\int_0^{\pi/2}\mathrm{round}\big(A\cos(x)\big)^2dx = \frac{\mathrm{round}(A)^2}{\pi/2}\,\mathrm{acos}\left(\frac{\mathrm{round}(A) - 0.5}{A}\right) + {\frac{1}{\pi/2}\sum_{k=1}^{\mathrm{round}(A) - 1}}k^2\Bigg(\mathrm{acos}\left(\frac{k - 0.5}{A}\right) - \mathrm{acos}\left(\frac{k + 0.5}{A}\right)\Bigg)\tag{6}$$

THD is calculated and minimized by the following Python script (continued), using Eqs. 4, 5, and 6:

def a_1(A):
    A = mp.mpf(A)
    return 2*(mp.floor(A + 0.5)*mp.sqrt(4*A**2 - (2*mp.floor(A + 0.5) - 1)**2) + mp.nsum(lambda k: k*(mp.sqrt(4*A**2 - (2*k - 1)**2)-mp.sqrt(4*A**2 - (2*k + 1)**2)), [1, mp.floor(A + 0.5) - 1]))/(mp.pi*A)

def MS(A):
    A = mp.mpf(A)
    return mp.floor(A + 0.5)**2*mp.acos((mp.floor(A + 0.5)-0.5)/A)/(mp.pi/2) + mp.nsum(lambda k: k**2*(mp.acos((k - 0.5)/A) - mp.acos((k + 0.5)/A)), [1, mp.floor(A + 0.5) - 1])/(mp.pi/2)

def STHD(A):  # Square of THD
    MS_1 = a_1(A)**2/2
    return MS(A)/MS_1 - 1

for m in range(2, 25):
    optimal_A = to_max_digits(lambda: mp.findroot(lambda A: mp.diff(STHD, A), approx_optimal_A(m)), 80, 100, double_digits)
    B = to_max_digits(lambda: a_1(mp.mpf(optimal_A)), 80, 100, double_digits)
    THD_optimal = to_max_digits(lambda: 10*mp.log10(STHD(mp.mpf(optimal_A))), 80, 100, dB_digits)
    THD_1 = to_max_digits(lambda: 10*mp.log10(STHD(2**(m-1)-1)), 80, 100, dB_digits)
    THD_2 = to_max_digits(lambda: 10*mp.log10(STHD(2**(m-1)-0.5)), 80, 100, dB_digits)
    print(str(m)+"&"+optimal_A+"&"+B+"&"+THD_optimal+"&"+THD_1+"&"+THD_2+"\\\\")

Figure 2. THD (dB) as function of unquantized amplitude $A$. The local maxima are at amplitudes that cause a narrow protrusion to poke out at the extrema of the quantized sinusoid.

Figure 3. Amplitude $a_1$ (blue solid line) of the fundamental frequency in the quantization of a cosine wave of amplitude $A$, with the identity line (orange dashed line) plotted for reference.

Table 2. Optimal $A$ by definition 2 for different $m\le24$ and the resulting THD, with the THD for some common choices of $A$ listed for comparison. Surprisingly, THD is minimized by the same $A$ that maximize SNR (Table 1), also when tested with much higher precision than what is shown here. If the signal is considered to be a sinusoid of amplitude $a_1$ instead of amplitude $A$, then SNR values are obtained by flipping the sign of the THD values.

$$\begin{array}{l|l|l|lll}m&\text{optimal }A&a_1&\rlap{\text{THD (dB)}}&&\\
&&A=\text{optimal}&A=\text{optimal}&2^{m-1}-1&2^{m-1}-0.5\\
\hline
2&1.26827949461530&1.17011951948679&-10.7629578&-10.1492078&-10.5728562\\
3&3.23800942121037&3.19552705145171&-18.7633376&-18.2533980&-18.3370094\\
4&7.21658597929407&7.19632525078776&-25.5045573&-25.1895549&-25.0267366\\
5&15.2007181331955&15.1907044658090&-31.8098475&-31.6159563&-31.3681507\\
6&31.1888756714257&31.1838597476996&-37.9357894&-37.8138122&-37.5642711\\
7&63.1800835394190&63.1775601103885&-43.9857175&-43.9071306&-43.6902709\\
8&127.173613625523&127.172343338577&-50.0049084&-49.9531877&-49.7783344\\
9&255.168894736361&255.168255766598&-56.0134634&-55.9788184&-55.8439127\\
10&511.165479188890&511.165158147298&-62.0199994&-61.9964656&-61.8952456\\
11&1023.16302205377&1023.16286093052&-68.0278654&-68.0117064&-67.9371468\\
12&2047.16126264484&2047.16118185697&-74.0380588&-74.0268736&-73.9726322\\
13&4095.16000722516&4095.15996674788&-80.0505958&-80.0428071&-80.0036830\\
14&8191.15911371601&8191.15909344706&-86.0651409&-86.0596937&-86.0316405\\
15&16383.1584789666&16383.1584688212&-92.0812822&-92.0774606&-92.0574286\\
16&32767.1580286428&32767.1580235662&-98.0986407&-98.0959534&-98.0816905\\
17&65535.1577094659&65535.1577069262&-104.116904&-104.115011&-104.104877\\
18&131071.157483397&131071.157482127&-110.135830&-110.134495&-110.127304\\
19&262143.157323352&262143.157322717&-116.155234&-116.154293&-116.149195\\
20&524287.157210089&524287.157209771&-122.174984&-122.174319&-122.170708\\
21&1048575.15712995&1048575.15712979&-128.194979&-128.194509&-128.191953\\
22&2097151.15707326&2097151.15707318&-134.215150&-134.21482&-134.213010\\
23&4194303.15703317&4194303.15703313&-140.235446&-140.235212&-140.233932\\
24&8388607.15700481&8388607.15700479&-146.25583&-146.25567&-146.25476\end{array}$$

The equivalence of the two optimality definitions holds even when numerical precision is increased significantly, here with $m=4$ at least to 200 decimal places, in Python (continued):

m = 4
mp.mp.dps = 200
mp.findroot(lambda A: mp.diff(RelMSE, A), 2**(m-1)-1+0.157)
mp.findroot(lambda A: mp.diff(STHD, A), 2**(m-1)-1+0.157)

which outputs numerically identical optimal values for $A$ for the two optimality definitions:

mpf('7.21658597929406951556806247230383254685067097032105786583650636819627678717747461433940963299310318715204551609940031954265317274195597248077934451075855527')
mpf('7.21658597929406951556806247230383254685067097032105786583650636819627678717747461433940963299310318715204551609940031954265317274195597248077934451075855527')

Limit $m\to\infty$

At large $m$ it becomes difficult to directly optimize $m$ numerically, so another approach is desirable. A Taylor approximation of the sinusoid about its peak, where it most differs from a linear function, is a quadratic polynomial. This can be used to analyze the effects of quantization in the limit $m\to\infty$ $\Rightarrow$ $A\to\infty$. The difference between the MS quantization error of a sinusoid with an amplitude $A$ and the MS quantization error $1/12$ of a linear function is proportional to (Fig. 4):

$$\begin{gather}\mathrm{MS} - \frac{1}{12} \propto f(a) = \int_0^{\sqrt{4a + 2}/2}\left((x^2 - a)^2 - \frac{1}{12}\right)dx + {\sum_{k=1}^{\infty}\int_{\sqrt{4a + 4k - 2}/2}^{\sqrt{4a + 4k + 2}/2}\left((x^2 - a - k)^2 - \frac{1}{12}\right)dx}
= \frac{1}{60}\left(\frac{\sqrt{4a + 2}}{16a^2 - 4a - 1} + \sum_{k=1}^\infty\left(\sqrt{4a + 4k + 2}\left(16a^2 + 4a\,(8k - 1) + 16k^2 - 4k - 1\right)\\
- \sqrt{4a + 4k - 2}\left(16a^2 + 4a\,(8k + 1) + 16k^2 + 4k - 1\right)\right)\right),\end{gather}\tag{7}$$

when the amplitude $A$ is an integer $\to\infty$ plus a real number $-0.5 < a \le 0.5$. The sum in Eq. 7 appears to converge, whereas leaving out the term $-\frac{1}{12}$ would result in the sequence of partial sums growing without bound, indicating that no matter which value of $a$ is chosen, $\left(\mathrm{MS} - \frac{1}{12}\right)/\mathrm{MS} \to 0$ and $\mathrm{MS} \to \frac{1}{12}$ as $m\to\infty$.

Figure 4. $f(a)$ and $f'(a)$ of Eq. 7. The shape of $f(a)$ looks identical to the shape of RRMSE for large $m$ in Fig. 1.

By symbolically differentiating $f(a)$ defined in Eq. 7 with respect to $a$ (Fig. 4) and by finding the zero of $f'(a)$ near $a = 0.157$, the optimal $a$ at $m\to\infty$ can be computed to quite a bit of precision, in Python (continued):

def f(a):
    a = mp.mpf(a)
    return (mp.sqrt(4*a + 2)*(16*a**2 - 4*a - 1) + mp.nsum(lambda k: (mp.sqrt(4*a + 4*k + 2)*(16*a**2 + 4*a*(8*k - 1) + 16*k**2 - 4*k - 1) - mp.sqrt(4*a + 4*k - 2)*(16*a**2 + 4*a*(8*k + 1) + 16*k**2 + 4*k - 1)), [1, mp.inf]))/60

def Df(a):  # Derivative of f(a)
    a = mp.mpf(a)
    return mp.sqrt(2)*(16*a**2 + 4*a - 1)/(12*mp.sqrt(2*a + 1)) + mp.nsum(lambda k: (mp.sqrt(4*a + 4*k - 2)*(16*a**2 + 4*a*(8*k + 1) + 16*k**2 + 4*k - 1) - mp.sqrt(4*a + 4*k + 2)*(16*a**2 + 4*a*(8*k - 1) + 16*k**2 - 4*k - 1))/(12*mp.sqrt(2*a + 2*k + 1)*mp.sqrt(2*a + 2*k - 1)), [1, mp.inf])

to_max_digits(lambda: mp.findroot(lambda a: Df(a), 0.157), 63, 83, double_digits)

which gives $a \approx 0.156936321399$. This can be used to create an approximation of the optimal $A$ as function of $m$, intended for large $m \ge 20$:

$A \approx 2^{m-1} - 1 + 0.156936321399 + e^{-1.2749819017
 - 0.3464088917 m},\tag{8}$

where the coefficients in the exponent were calculated by a linear fit at $m \in \{21, 22\}$ to a linearization of optimal $A$ values of Table 1 or equivalently Table 2. The error from using the approximation is shown in Fig. 5.

Figure 5: Absolute error in the approximation of optimal $A$ by Eq. 8 as function of $m$. For $m \in \{21, 22\}$, which were used for fitting, and for $m \in \{23, 24\}$, the absolute approximation error was less than $10^{-8}$.

Out of curiosity, I also computed the worst-case $a$ that gives the largest MS quantization error as $m\to\infty$, by finding the zero of $f'(a)$ near $a = -0.43$, which turned out to be at $a \approx -0.433510875868$.

Conclusion

The two definitions of optimality appear to be equivalent, to convincing numerical precision. As $m\to\infty$, the optimal value of $A$ approaches approximately $2^{m-1}-1+0.156936321399$ (or more accurately for large finite $m$ the approximation of Eq. 8) and the quantization error reduction (by def. 1 SNR or def. 2 THD) in dB from choosing the optimal value approaches zero, compared to choosing another large value such as $A = 2^{m-1}-1$ or the nearly worst-case choice $A = 2^{m-1}-0.5$.

At the optimal amplitude $A$, the least squares (LS) sinusoid is of the same frequency and phase as the sinusoid being quantized, but has a somewhat lower amplitude $a_1$ given in Table 2. This is somewhat counterintuitive. To minimize THD (or to maximize SNR with the sinusoid being approximated as the signal) of approximating a sinusoid of amplitude $a_1$ using a waveform of $m$-bit numbers with values in range $-2^{m-1} \ldots 2^{m-1}-1$ that is constructed by quantizing (by rounding to the nearest integer) a sinusoid of amplitude $A$, one must choose optimal $a_1$ and $A$ that are not equal.

The results are applicable to continuous time without sampling, or to sampling in the limit $f\to f_s$, where $f$ is the sinusoid frequency and $f_s$ is the sampling frequency. In general, optimality is not preserved by sampling. Random sampling with random times of the samples will preserve definition 1 SNR optimality, if the distribution of the sinusoid phase (modulo $2\pi$) at the samples is uniform. Also, for irrational $f/f_s$, the same optimality is preserved for the same reason, see equidistribution theorem. Definition 1 SNR optimality vanishes with rational $f/f_s$, because the distribution of sinusoid phases will not be uniform.


Can Normalized Cut be used to find connected components in binary images? I'm unsure whether Normalized Cut can or can't be used to find connected components in binary images. I would argue that it could, but could not validate this argument.

Lamar, first welcome to StackExchange. Laplacian of a graph has very strong links to the topological properties such as connected components. For simple scenarios, the graph Laplacian, $L=I-D^{\frac{-1}{2}}AD^{\frac{-1}{2}}$ constructed using the degree matrix $D$ and the adjacency matrix $A$, can be processed to result in connected components. First, the number of connected components in the graph is the dimension of the nullspace of the Laplacian and the algebraic multiplicity of the $0$-eigenvalue. Then, for a graph with multiple connected components, $L$ is a block diagonal matrix, where each block is the respective Laplacian matrix for each component, possibly after reordering the vertices.

Simply look at the eigenvalues of the Laplacian matrix and identify clusters (e.g. Spectral clustering). This will give you the individual components if you form the graph with correct adjacency matrix and degrees.

I'm not sure, but if you directly like to use normalized cuts for finding the CC: Simply form the graph ($1$s are adjacent and $0$s are adjacent in themselves). Run the NC and select the segments which do not belong to the background. These should be the connected components.


How to set covariance for Bivariate Logistic Distribution This is the logistic distribution of single random variable (taken from Wikipedia).
x = random variable
mu = mean of all random variables
s = variance.

Now, I want to do a Bivariate logistic distribution (having two random variables x1 and x2). My dataset is going to be image pixel values!
When I find the covariance of two random variables, it turns out to be a 2 x 2 square matrix but I need a single number! How to actually compute a bivariate model?

The covariance of any multivariate distribution is always a matrix. Maybe it is of help to you to know what this matrix contains:

-In general, each entry i,j of the matrix gives you the correlation between the RVs i and j.

-So a diagonal element i,i will simply be the variance of the RV i.

-If RV i and j are independent, then the non-diagonal entry i,j is zero (they are uncorrelated). The reverse of this statement is not generally true, but it is for Gaussians iirc.

I don't know what you are trying to do, but perhaps the trace of the matrix or its determinant is of use to you?


Estimation of CCD/CMOS sensor noise I would like to know what is the basic procedure to estimate the noise level of a CCD/CMOS sensor from an image/sequence of images acquired. I am going to use MATLAB to perform this operation.

What I did is I acquired a certain number of images with the shutter on, so the images are basically black images. The cameras I tested so far are all colour cameras, so there's a Bayer filter in front of the sensors. I have three channels R,G,B in the image output.

I was thinking that I can compute the average pixel value and the two-dimensional RMS to have an estimation of the noise. But I don't know if this is the recommended approach.

In addition, I guess that the noise estimation with the shutter on does not relate much with the noise level with the shutter on.

What do you think, do you think that the lens aperture also plays a significan role?

Instead of capturing with shutter on, I suggest you to capture a picture in low-light and set your exposure time as minimum as possible and the reduce the aperture. Capture a neutral surface(where R,G,B are equal,preferably neutral gray) so when averaging the entire frame you will know noise level change by particular color.Do your calculations on Bayer RAW instead of jpeg because RAW is directly from sensor without any extra processing.

  Do you think lens aperture will play a significant role?

Yes.Lens aperture will play a significant role. If you are reducing the aperture it means the amount of light reaching the sensor is less and it has to be compensated by either increasing exposure time or the gain(Increasing the gain increases the noise).


Do we assume modulating and carrier signal uncorrelated in AM modulation? I am new to communication and studying amplitude modulation. Let us assume an amplitude modulated wave given by
$$\big(m(t)+ A\big) \cos(2\pi f_c t)$$
Now we have formula for efficiency as
$$ \eta=\text{useful power / Total power}=\frac{\overline{m(t)^2}/2}{A^2/2+\overline{m(t)^2}/2}$$
I am not able to understand why total power is not equal to $$\overline{(m(t)+ A)^2}$$ where $\bar{.}$ denotes average power. I mean where does the term for $A^2\overline{m(t)^2}$ go in the denominator of the efficiency expression. Do we assume modulating signals and carrier signal to be uncorrelated?

The total power of an amplitude modulated signal is
$$\begin{align}\overline{s^2_{AM}(t)}&=\overline{\big(A+m(t)\big)^2\cos^2(2\pi f_ct)}\\&=\frac12\overline{\big(A+m(t)\big)^2}+\frac12\overline{\big(A+m(t)\big)^2\cos(4\pi f_ct)}\tag{1}\end{align}$$
The second term on the right-hand side of $(1)$ is zero if $m(t)$ is a lowpass signal, and if $f_c$ is sufficiently large (i.e., greater than the highest frequency in $m(t)$).
Consequently, we have
$$\begin{align}\overline{s^2_{AM}(t)}&=\frac12\overline{\big(A+m(t)\big)^2}\\&=\frac12\big(A^2+2A\cdot\overline{m(t)}+\overline{m^2(t)}\big)\tag{2}\end{align}$$
The term $\overline{m(t)}$ is the DC value of the message signal, which is assumed to be zero. With this assumption we finally arrive at
$$\overline{s^2_{AM}(t)}=\frac12\big(A^2+\overline{m^2(t)}\big)\tag{3}$$
In sum, the assumptions for which $(3)$ is valid are that $m(t)$ is a lowpass signal, that the carrier frequency is sufficiently high, and that $m(t)$ has a DC value equal to zero. It is usually straightforward to satisfy all of these assumptions to a very high degree of accuracy.


Sub pixel shift causes weird results I am trying to shift an image on a subpixel level, with the method described here:
http://www.mathworks.com/matlabcentral/fileexchange/18401-efficient-subpixel-image-registration-by-cross-correlation/content/html/efficient_subpixel_registration.html

The results are however disappointing. Since I am actually programming this algorithm for the rotation of images, where I resample inbetween pixels, this problem keeps on bugging me: http://spin.ecn.purdue.edu/fmri/PDFLibrary/TongR_MRM_1999_41_253_256.pdf

Can somebody tell me how these strange artifacts can be resolved? According to multiple scientific articles I found, it should be possible to resample on other grids with satisfactory results. My minimal code example (in python) is below, along with the resulting images, which are not satisfactory at all.

Does someone know what causes this effect, and how it can be resolved to reproduce the results others have?

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

N = 2048

img=mpimg.imread('koe.png') #http://www.clipartbest.com/cliparts/RiA/yXL/RiAyXLriL.png

f = np.zeros((N,N))
f[:img.shape[0],:img.shape[1]] = img[:,:,3]

def fft(x):
    return np.fft.fft2(x)

def ifft(x):
    return np.fft.ifft2(x)

exp = np.exp
pi = np.pi
j = 1j

def s(r):
    plt.figure()
    #plt.imshow(np.real(r), vmin=0.0, vmax=1.0, interpolation='none')
    plt.imshow(np.absolute(r), interpolation='none', vmin=0.0, vmax=1.0)
    plt.colorbar()

x = np.tile(np.arange(N),(N,1)).T
y = x.T
result = ifft( fft(f)* np.exp(1j*2*np.pi*(-0.5*x/N-0.5*y/N)) )
s(result)
s(f)
plt.show()

The solution lies in the use of the x's. You need to convert the x's you use to lower frequencies, which can be either positive or negative. So only use x's between -N/2 and N/2 and the problem is solved.

def convertfreqs(x):
    return np.mod((x + N/2),N) - N/2

This works because it is the higher frequencies which cause the signal to change heavily inbetween pixels. Higher frequencies means higher x's.


Stabilizing the inverse transform of a system Currently I am trying to understand an unknown mechatronic system. I used white noise as an input signal and measured the output of the system. Then I estimated a second order SISO transfer function (Output/Input) in the frequency domain using Matlab.
sys=Out/In
h = idfrd(sys, freq*2*pi, Ts);
np = 2;
nz = 2;
Gs = tfest(h, np, nz)
opts = c2dOptions('Method', 'tustin', 'PrewarpFrequency', 2*pi*120);
Gz = c2d(Gs, 1/4000, opts);

$$
Gs=\frac{3.765\cdot s^2 - 528.7 s - 5.113 \cdot 10^4}{s^2+183 \cdot s + 4.117 \cdot 10^5}
$$
A bode Plot of the derived second order system $Gs$ (blue) and the measured data (orange) shows good agreement in amplitude and phase above $60\;Hz$.

Further, a pole-zero plot of $Gz$ shows that all poles are stable (inside unit circle of z-plane), but that there is one zero outside of the unit circle of the z-plane.

Now, I want to counteract the behaviour of the mechatronic system in a way that I achieve a unity output. So by modifying the input to the system I want to basically get rid of the system characteristics and achieve a pure unit transfer behaviour.
IGz=1/Gz; % so that IGz * Gz = 1

Mathematically, I need to multiply the input signal with the inverse transform of the system. Here a bode plot of the system $Gs$ (blue) and its inverse transform $IGz$ (orange) basically shows what I want to do.

But because of the zero outside the unit circle, the inverse transfer function now has a pole outside the unit circle and thus is unstable. Therefore I wonder how I can counteract the behaviour of the system while not getting an unstable system?

Attempt
So far, I have tried to replace the unstable pole by reflecting the pole inside the unit circle as advised here: $(z-p)$ by $(z-1/p)$ (as the unstable pole is purely real I don't need to do the complex conjugate). This keeps the amplitude response, but changes the phase response (blue: unstable $IGz$, orange: stabilized $IGz$):

However, I would like to keep the phase response as close as possible to my inverse transfer function as well. But, as I don't care so much about what happens below $60\; Hz$ I wonder if there is an approximation of my inverse system which is good above $~60\;Hz$ in amplitude and in phase, but can deviate below $~60\;Hz$ and is stable. I wonder if there is an established way how one would derive such a system?
Here is what I have tried so far.
pIG=pole(IGz);
zIG=zero(IGz);

z=tf('z',1/4000);
Test=((z-zIG(1))*(z-zIG(2)))/((z-1/pIG(1))*(z-pIG(2)));
Fac=evalfr(IGz,100)/evalfr(Test,100);
Test1=Fac*Test;
bode(IGz,Test1)

The result looks good at higher frequencies, however the phase could be better below $100\;Hz$ or so (blue: unstable $IGz$, orange: stabilized $IGz$ with phase kept roughly the same above $100\;Hz$):

Edit
I just looked at the impulse responses. The system $Gs$ of course is oscillating, but decaying.

However, the impulse response of the inverse system is, of course, unstable.

The impulse Response of the second stabilized system is still overshooting with high amplitude. I wonder if this works, since I, of course, am limited in output w.r.t. the controller.

With the input $u$ to the controller and the output $y$, I want to do the following $u \rightarrow IGz (controller) \rightarrow y \rightarrow Gz (mechatronic\; system) \rightarrow y_2$. During all this, the output $y$ of the controller can't get too big and it should obey $\frac{y_2}{u}=1$.
How can I best come up with an $IGz$ that achieves this?

Like you mentionned, you cannot cancel a right-half-plane zero (or a zero outside the unit circle) by placing a pole on it. A unstable pole in your compensator will make the command of your controller unbounded (i.e. it will reach infinity).

There are no ways to cancel a right-half-plane zero. It's sometimes possible to "remove" a right-half-plane zero by chansing a sensor position.

From Murray's book http://www.cds.caltech.edu/~murray/books/AM08/pdf/fbs-limits_18Aug2019.pdf

  The existence of the right half-plane zero can be removed if we choose to measure
  the location of the vehicle by the position of the center of the rear wheels instead
  of the center of mass.

Edit : A right-half plane zero can also be caused by a time delay. Delays to computation can be reduced, but some other delays cannot be avoided.

Edit 2 : I will try to answer your other questions later.


Why having a ridge like Ambiguity Function implies lack of Doppler resolution? I was wondering why if the AF (Ambiguity Function) has a ridge like form it has low Doppler resolution.

If you examine a contour plot of the Ambiguity function - the lines connect equal values of the Ambiguity function. This means that if you travel along one of these lines (i.e. keeping the same value of the function) then your system will not be able to distinguish between the signal that have this combination of Doppler and time delay. Note that the values for time/Doppler resolution depend on what value of the Ambiguity function you choose i.e. which contour line you decide on. I'm saying just pick one for the purposes of illustration.

If you look the Ambiguity function for a CW pulse of 2 secs duration:

You see that if you choose the lowest contour value, that we would be unable to distinguish between time delay of -1.5 secs and 1.5 sec, and our Doppler resolution ranges from -0.5 Hz to 0.5 Hz. If you were use a pulse of longer duration the Ambiguity function would be wider (left to right) and shorter (top to bottom). This means your time resolution would be worse, but your Doppler resolution would be better. The opposite is true if you used a shorter pulse.

If we examine the Ambiguity function for a LFM chirp pulse:

Now we see that the ambiguity exists for a combination of Dopplers and time delays (chose one of the contours that is essentially an ellipse).Previously, for the CW pulse, across a range of time delays the ambiguity was mostly constant for the same Doppler value. The advantage of the LFM pulse over a CW pulse, is that I can use the same pulse length but achieve better time resolution. You could also use a longer pulse and achieve better time resolution. The longer pulse means you can get more energy on the target and improve the signal to noise ratio.


how to detect corner in a digital image using MATLAB? I am designing an algorithm for image processing. I am using MATLAB for the processing of the images(video). How I can detect the position of the corners in the image?

In the case of video, the best approach to use is the ShiTomasi Corner Detector. (detectMinEigenFeatures). Now, if you want to detect corners that contains some invariant featrues, you should use SURF, BRISK, ORG, AKAZE, FREAK, etc...

Example:

im = rgb2gray(imread('image.jpg'));
[r c] = corner_ST(im,20);


Mutual information of $ \infty $-PAM I'm trying to compute using matlab the mutual information for an $ \infty $-PAM input (the limit of a very dense PAM constellation) for a range of snr and I got stuck.

I'm working with a real-valued scalar Gaussian channel with input power constrained to unity and standard Gaussian noise $ N \sim (0, 1) $. The output of the channel is given by,
$$ y = \sqrt{ \gamma } x + n $$

I know that under the power input constraint the input is uniformly distributed on $ \left[ -\sqrt{ 3 }, \sqrt{ 3 } \right] $.

I found that in HIGH-SNR regime, there exists a gap between the Shannon capacity and the capacity achieved with $ \infty $-PAM around 0.254 bits. I know this is called the shaping loss and represents the loss of using an uniform input rather than a Gaussian input.

I know how to compute that. But, I got stuck trying to compute the capacity for a range of snr or $ \gamma $. I tried using the formula,
$$ I(x;y) = h(y) - h(n) $$
with $ f_{y}(y) = \int f_{x}(x)f_{y|x}(y|x) dx $ and I got
$$ f_{y} (y) = \frac{ 1 }{ 2 \sqrt{ 3 \gamma } } \frac{ 1 }{ 2 } \int_{ \frac{ y - \sqrt{3} }{ \sqrt{2} } }^{ \frac{ y + \sqrt{3} }{ \sqrt{2} } } \frac{ 2 }{ \sqrt{ \pi } } e^{ -t^{2} } dt \\
= \frac{ 1 }{ 2 \sqrt{ 3 \gamma } } \frac{ 1 }{ 2 } \left[ \text{erf} \left( \frac{ y + \sqrt{ 3 \gamma } }{ \sqrt{ 2 } } \right) - \text{erf} \left( \frac{ y - \sqrt{ 3 \gamma } }{ \sqrt{ 2 } } \right) \right] $$

I tried to compute $ h(y) = \int_{-\infty}^{\infty} f_{y} (y) \log_{2} f_{y} (y) $ and the integral does not converge. I tried with int$(\cdot)$ and integral$(\cdot)$ methods in matlab.

I also tried approximating the $ \text{erf}(\cdot) $ function with Taylor series expansion
$$ erf(x) = \frac{ 2 }{ \sqrt{ \pi } } \sum_{ n = 0 }^{ \infty } \frac{ (-1)^{n} x^{ 2n + 1 } }{ n!(2n + 1) } $$
but it didn't work.

How can I obtain a more manageable integral? Is my approach correct?

I appreciate any help!! Thanks.

I got an answer here,

https://math.stackexchange.com/questions/1554659/mutual-information-for-a-continuous-uniform-distribution

The $ \text{erf}(x) $ function approximates to 1 with an error less than 0.0001 for $ x \geq 3 $. I plotted $ f_{y} (y) $ and I saw that the most of probability is within a band denpending on $ \sqrt{ 3 \gamma } $. Moreover, taking the limit $ \lim_{ y\to\infty } g(y,\gamma) \log g(y,\gamma) = 0 $. Then, I computed

$$ h(y) = \frac{1}{2} \log \left( 24 \gamma \right) - \frac{ 1 }{ 2 \sqrt{ 3\gamma} } \frac{1}{2} \int_{ - \sqrt{ 3\gamma } - a }^{ \sqrt{ 3\gamma } + a } g(y,\gamma) \log_{2}( g(y,\gamma) ) dy $$

with

$$ g (y,\gamma) = \text{erf} \left( \frac{ y + \sqrt{3\gamma} }{ \sqrt{2} } \right) - \text{erf} \left( \frac{ y - \sqrt{3\gamma} }{ \sqrt{2} } \right) $$

using $ quadgk(\cdot) $ function in matlab.


Instantaneous Frequency Estimation by Hilbert Transform - Theoretical Justification and Proof I would like to better understand why the instantaneous frequency estimation by Hilbert transformation works (and especially why it doesn't work / lead to precise results in many cases).

The motivation is to estimate signal $x(t)$ by decomposing it into an amplitude envelope $m(t)$ and phase of cosine $\omega_c (t)$ (or carrier waveform):

$$x(t) = m(t)  \cos\left(\omega_c (t)\right)$$

Now, assume that $x(t)$ indeed is a result from such a process.

Questions:

1) There are two "parameters" to be estimated for any $t$, as such some constraints are needed. What are the constraints regarding $m(t)$ and $\omega_c (t)$ that are selected when applying the Hilbert transform decomposition?

2) Is there a proof available somewhere that given the constraints, the estimation indeed finds the correct amplitude envelope and carrier (for continuous and also discrete case)?

See my earlier comments here: Meaning of Hilbert transform

Common fractal noise isn't an analytic signal (infinitely differentiable).  And a Hilbert transform re-creates the imaginary component of an analytic signal if you have the real component of the complex analytic signal (which one rarely has from real-world data).

However, sufficiently band-pass filtered data might be similar to a finite length segment of something that looks like an infinite length analytic signal (e.g. is from a source whose behavior can be modeled or estimated by a 2nd order linear differential equation).


Butterworth band-pass filter transfer function I have some complex data in the frequency domain and want to multiply by an appropriate function to get complex data in the frequency domain that is appopriately filtered.

I think such a function should be, in terms of frequency:

$$ 1 - \frac{1}{\sqrt{1 + (\frac{fW}{f^2 - f_0 ^2})^{2n}}} $$

where W is the width of the frequency band to pass, and $f_0$ is the center frequency.

Is the above correct, or what would be the correct function to accomplish this?

okay, the magnitude frequency response for a general $N$th-order Butterworth lowpass filter is:

$$ \Big|H_\mathrm{LP}(f)\Big|^2 = \frac{1}{1 + \left(\frac{f^2}{W^2} \right)^N} $$

for bandpass filtering, centered at frequency $f_0$, this substitution is made to the LPF:

$$ f \leftarrow \frac{f}{f_0} - \frac{f_0}{f} $$

so then it comes out as:

$$ \Big|H_\mathrm{BP}(f)\Big|^2 = \frac{1}{1 + \left(\frac{\left(\frac{f}{f_0} - \frac{f_0}{f}\right)^2}{W^2} \right)^N} $$

blast that out and you'll have the frequency response for an $2N$th-order Butterworth bandpass filter.


Dealing with spectral leakage between two frequency bins I'm looking for two signals embedded in very noisy data. I know they are close to periodic and integral with the run, but are rather close (30s and 21s during a 210s measurement), and our sampling rate is low (1.25s). We have multiple observations to work with per sample.

I am wondering how I might deal with the problem of spectral leakage: it is evident in the FFT of our samples. Since our two signals of interest are so close together (10 and 7 cycles per measurement), it may well be that spectral power is leaking from the stronger of the two signals into the weaker one.

Will zero-padding with a window (rectangular or hamming) allow us to define a boundary between the two signal? Would another approach would be to concatenate these observations, driving the periodicities further apart (20 and 14)?

Thanks for your help, joseph

Yes, concatenating the data will greatly reduce the spectral leakage problem.  There will still be spectral leakage of course, but there will be more bins in between your signals-of-interest, so they won't leak as much into each other.

And yes, windowing should help too.  You are already using a "rectangular" window.  I would suggest trying other windows such as Hamming.


How do I model a system to use with a Kalman Filter? So I have an accelerometer, gyroscope, magnetometer and GPS. I would like to use a Kalman filter to optimally measure speed, position, acceleration and orientation. I have done research and I need to model the system mathematically, however this is the part I have got stuck on.

Can anyone give me some pointers on how to proceed here, do I need to try and model the system using differential equations and then derive the parameters for my $\mathbf A, \mathbf B$, and $\mathbf H$ matrices from there?

Any guidance would be great, thanks.

The following article explains the questions you have on modeling the system mathematically, in a very lucid way with examples and equations.

http://www.anuncommonlab.com/articles/how-kalman-filters-work/index.html


Calculating the bit period I want to calculate the bit period in my RF packets that I transmit them with 50Kbaud (Symbol rate=bit rate in my case) , how can I calculate the bit period?

50 KBaud means the symbol rate is 50 KHz (and in the OP's case the bit rate with 1 bit per symbol). The units of Hz is 1/sec. So the duration of each bit (the bit period) is simply

$$t_p = \frac{1}{50,000} = 20 \mu\text{s}$$


White noise vs. delta pulse and Ultraviolet catastrofe Everybody explains that white noise has all frequencies equally strong. But, this immediately means that

Ultraviolet catastrofe inevitably happens if power > 0 stays constant at any frequency and, what is similarly unacceptable,
White noise is identical to single Dirac impulse since delta pulse is a constant in the Fourier basis). Note that constant is the opposite to the notion of noise.

Is average over time is implied, like you obtain uniform distribution, or what?

Why answers are not satisfactory

Answers say that white noise does not mean that signal spectrum is limited. It means that that power spectrum is flat, which also means that there is no correlation between samples so that it is correlation matrix, which is delta function (and I do not know what is its fourier transform is, which has the infinite spectrum).

I see a several problems with that. At first, the power spectrum is simple a square of the spectrum. If you say that you prevent the Utraviolet catastrophe by cutting all frequencies above some threshold $f_{max},$ you cannot have a flat spectrum anymore.

Secondly, I understand that you can have a mean and variance of a uniform distribution which has value $v$ in $[a,b]$ and 0 outside. But what is a mean and variance of a perfectly flat power spectrum? Ok, mean might be zero if you admit negative frequencies. But you say that variance is $\sigma^2$. How is that?

Lastly, I have determined that rapid changes are less likely in the the limited spectrum signal, which means that samples tend to each other, which means that they are correlated. Ok, might be you say that they are correlated but pink nose does not say if they are positively or negatively correlated, so samples are not correlated in case of pink noise. Ok, this is great. But we have just concluded that pink noise is white (or can be white). Is it right?

I also see Wikipedia saying that white noise can be Gaussian, which means that the samples are normally distributed. This means that they will tend to each other, like in the pink noise.

As Jason was saying in the comments, the power spectral density of white noise is flat.  This is equivalent to saying that the autocorrelation of white noise is a delta dirac function (i.e. that there is no correlation, positive or negative, between one noise sample and another), not that the noise itself is a delta dirac function.


inquiries for writing continuous wavelet transform codes manually I want to write continuous wavelet transform  codes manually by matlab. And I want to use complex morlet function. Here are some background:
Continuous wavelet transform definition: $C(S,T;f(t),\psi (t))=\frac{1}{\sqrt{S}}\int_{S}^{b}f(t)\psi^{\ast }(\frac{t-T}{S})$
S is scale vector.for example 1:60;$T$ is time slid.
And $\psi(t)= \frac{1}{\sqrt{\pi f_{b}}}e^{2\pi f_{c} t}e^{\frac{-t^{2}}{f_{b}}}$

psi = ((pi*fb)^(-0.5)).*exp(2*1i*pi*fc.*t).*exp(-t.^2/fb); % for example fb=15;fc=1;

my discrete signal has $N$ points. discrete version of that integral for the first element in the scale vector  is:  for s(i)
$C_{T}(s)=\sum_{n=0}^{N-1}f(n)\psi^{*}(\frac{T-n}{S})$ this must calculate for all scales,1:60;
at the end it must return me a complex matrix with $N\times S$ size.
I am confused about wrting this codes. Anyone can help?
P.S I don't want to uses matlab function conv2 to calculate that convolution.
Thanks in advance.
 here is my first attempt, but it's not working at all.

%% user CWT
clear all
N=300;                %sample point numbers
t=linspace(0,30,N);
%% signal
x=5*sin(2*pi*0.5*t);  % signal with freq of 0.5 HZ
%% cwt
    fc=1;fb=15;
% psi=((pi*fb)^(-0.5)).*exp(2*1i*pi*fc.*...
%     t).*exp(-t.^2/fb);
%% convolution Psi([N-n]/S)*x(n) so we calculate convolution(psi(n/s),x(n))
    for s=1:60   %scale vector s=[1:1:60]
for i = 1:N      % number of discrete times
       for k = 1:i
               if ((i-k+1)<N+1) && (k <N+1)
                PSI = ((pi*fb)^(-0.5)).*exp(2*1i*pi*fc.*...
       (t/s)).*exp(-(t/s).^2/fb);
                        c(i,s) = c(i)+ x(k)*PSI(i-k+1);
                end
       end
end
    end

Here is the Python code that I use for CWT. You'll need a short 16bit mono 44.1Khz test.wav soundfile as input.

In order to compute the CWT, we need to compute the convolution between the input x[n] and the morlet wavelet. An efficient way to do this is to do ifft(x_ft * morlet_ft)
( I understand this as $f * g(t) =\frac{1}{2 \pi} \int_{-\infty}^{+\infty} \hat{f}(w) \hat{g}(w) e^{i t w} d w$ : it is more efficient to do a multiplication of the fourier transforms and THEN to do an inverse fourier transform, than to compute a convolution directly)

from scipy.io.wavfile import read
from pylab import *
import matplotlib.pyplot as plt
import numpy as np
import scipy

# read the file
(sr, samples) = read('test.wav')
x = (np.float32(samples)/((2 ** 15)-1))
N = len(x)
x_fft = np.fft.fft(x)

# scales
J = 200
scales = np.asarray([2**(i * 0.1) for i in range(J)])

# pre-compute the Fourier transform of the Morlet wavelet
morletft = np.zeros((J, N))
for i in range(J):
  morletft[i][:N/2] = sqrt(2 * pi * scales[i]) * exp(-(scales[i] * 2 * pi * scipy.array(range(N/2)) / N - 2)**2 / 2.0)

# compute the CWT
X = empty((J, N), dtype=complex128)
for i in range(J):
  X[i] = np.fft.ifft(x_fft * morletft[i])

# plot
plt.imshow(abs(X[:,scipy.arange(0,N,100)]), interpolation='none', aspect='auto')
plt.show()

Here is the scalogram plot that this code gives :


Filter design to realize Cauchy product I come from Computer Science so please pardon for my possibly wrong terminology.

I need to design a filter which has coefficients
$$h_0, h_1, \ldots, h_n, \ldots \quad\text{such that}\quad h_0 > h_1> \ldots > h_n > \ldots$$

The input is $x_n$. The output is $y_n$. The relationship between $h_n$, $x_n$ and $y_n$ is a convolution:
$$y_n = h_0 x_n + h_1 x_{n-1} + \ldots + h_n x_0$$

The problem is that I have to carry out a convolution to compute the output $y_n$ upon the arrival of every new input $x_n$. So I have to keep all $x_1, x_2, \ldots, x_n, \ldots$ in memory and do $n+1$ multiplications and $n$ additions every time.

I would like to ask whether there is any way to design the coefficient $h_n$ such that I can compute $y_n$ more efficiently. I know of one solution, which is $h_n=\lambda^n$ with $0<\lambda<1$ but this function decays too fast.

Any idea to approximate this procedure is also welcomed.

You could implement some convolutions as recursive filters. Consider moving average as an example, to calculate the output signal, instead of summing over the averaging window, you could add the new input and subtract the part of input which exited the window from the output signal to obtain the output of next time step. Also some filters could be implemented as lattice.

But there is not a general method to implement all filters recursively or as lattice.


Design of FIR and IIR filters in Matlab I need to create 4 different pass-band filters. What remains the same throughout is the following:

Sampling Frequency = $3000Hz$

Pass-band Frequency = $900Hz-1200Hz$

Now, these are the filters I am supposed to create:

$24^{th}$ order FIR using windowing method with Hamming window.
$24^{th}$ order FIR  using frequency sampling.
$8^{th}$ order IIR using bilinear transform with Butterworth prototype analogue filter
$8^{th}$ order IIR using impulse invariance with Butterworth prototype analogue filter.

My Attempts:

Normalised frequencies are $0.6$ and $0.8$ and the filter is 24th order  hence

b = fir1(24,[0.6 0.8]);

The Desired frequency response is
$$H(\Omega)=\exp^{-12j\Omega}$$ for $$\frac{3\pi}{5}\leq\Omega<\frac{4\pi}{5}$$
and
$$\frac{6\pi}{5}<\Omega\leq\frac{7\pi}{5}$$
With $$\Omega=\frac{2\pi k}{25}$$
Hence I have that $H(\Omega)\neq0$ when $k=8,9,16,17$

The function in Matlab that creates this filter is

B=fir2(N,f,m)

Where $f$ is the vector with sampled frequencies and $m$ is the magnitude at those frequencies. But for my case I have frequencies higher than $1$ at it says that $f$ starts with $0$ and ends with $1$. I kind of guessed the outcome and produced this:

f=[0 0.64 0.72 1];
m=[1 exp(-12*j*0.64) exp(-12*j*0.72) exp(-12*j)];
b=fir2(24,f,m)

I have a feeling that this is not correct.

For the bilinear IIR I went to matlab's "butter" function.

[b,a] = butter(4,[0.6 0.8],'stop');

For this part I know I need to use both 'butter' and 'impinvar' together but I don't know how.

EDIT / ATTEMPT:

[b1,a1] = butter(4,[2*pi*900 2*pi*1200],'s');
[bz,az] = impinvar(b1,a1,3000);

Could someone tell me if my attempts 1-3 are correct and help me out a little with 4? Thank You.

Your call to fir1 looks OK, but you should check the result yourself. You could probably answer several of your other questions by carefully reading the mathworks documentation of the respective functions. E.g., for fir2 you can read that

  f is a vector of frequency points in the range from 0 to 1, where 1 corresponds to the Nyquist frequency. The first point of f must be 0 and the last point must be 1. f must be sorted in increasing order. Duplicate frequency points are allowed and are treated as steps in the frequency response.

and

  m is a vector containing the desired magnitude response at each of the points specified in f.

So you shouldn't define m as a complex vector; the routine knows the phase lag because it's just half of the filter order. The vector f should be defined as follows:

f = [0 .6 .6 .8 .8 1]

Since this is homework, I leave the definition of the magnitude vector up to you.

Your call to butter looks OK, apart from the fact that you used stop as the filter type. You want a band pass filter, not a band stop filter, so just leave out the additional filter type argument, and the routine will by default design a band pass filter (because there are two edge frequencies in the frequency vector).

For the impulse invariance method, you first have to design an analog Butterworth filter (using butter with last input argument s for analog filter), and then use impinvar to create the corresponding digital filter. The mathworks documentation is very clear, and I think you will succeed after having studied it. Otherwise, add to your original question a very specific question where you're stuck.


Working out integration of basic signal functions As the title suggests , I'm having problems working out the integration of some basic signal functions. Shown below are the questions and my workings. For some reason I'm just unable to get the answer.

For a) and c) , I know that for u(t) is basically 0 before t=0 , thus I changed the limits for the integration. I did the same thing for u(t-1)

Thanks for the help!

So you've already figured out the $u(t-T)$ trick:  the unit step (delayed) just changes the limits of integration.

The other trick is that $\delta(t-T)$ just "samples" the value of the rest of the integral's argument. So, for example, $\int x(t) \delta(t-T) dt = x(T)$ where $T$ is somewhere within the integral's limits.


Discrete Fourier style transform with non-sinusoidal/arbitrary waveform components A discrete Fourier transform will produce frequency samples sufficient to recreate the original time sampled signal from sinusoidal waveforms.
That works great as long as the original signal is actually composed of sinusoidal waveforms. Let's assume that we know a signal is composed of non-sinusoidal waveforms, under which assumption its decomposition may be simpler, such as triangular waveforms. Does a discrete Fourier style transform exist which will decompose the sampled signal into triangular, rectangular or other waveforms? Does such a transform exist for any arbitrary waveform?
For example, here is the FFT response of a 150 Hz pure sin wave sampled at 1024 Hz (courtesy numpy and matplotlib):

Here is the FFT of a 150 Hz pure sawtooth, again sampled at 1024 Hz:

It's desirable here to have a similarly "clean" transform using the assumption that the signal is composed of sawtooth waveforms. So my question, does a transform exist that would decompose the signal into some non-sinusoidal waveforms?

Any complete orthogonal set of functions can be used to represent any discrete function—just like trigonometric functions can be used to represent any continuous function in Fourier analysis. You can use "triangular, rectangular or other waveforms", only mind those represent an orthogonal set. There is also a "wavelet transform", see Wikipedia. For transforms using Walsh functions, the fast implementations exist: Fast Walsh-Hadamard transform.
For visualization of the concept, see images in Wolfram MathWorld pages:

downloaded from https://mathworld.wolfram.com/WalshFunction.html


Non-causal FIR filter in the feedback loop I have a feedback loop with a transfer function $H(z) = \sum_{i=0}^{L-1} h(i) z^{-i}$. Is there a way to make this FIR filter non-causal? If it was a feed-forward filter, we could simply do so by buffering data. But can it be done in a feedback loop?

Is there a way to make this filter non-causal?

Remember that non-causal filters aren't possible to interpret in any case, because "non-causal" literally means there's output caused by input that comes later.
Just because you can buffer something if latency doesn't matter doesn't mean you've actually built a non-causal system - your buffering adds delay back in that actually has to be at least as long as the anticausal part of your system.
So, no. This makes no sense. You'll have to make your H(z) causal (which isn't hard at all), and deal with the resulting group delay (which might be problematic).


Does $ Ev(\cos(4\pi t)u(t)) $ signal has period? According to Oppenheim's signal and systems solution book and other solutions like this

$$
Ev(\cos(4\pi t)u(t)) = 1/2 \times ( \cos(4\pi t)u(t) + \cos(4\pi t)u(-t)) = 1/2(\cos(4\pi t))
$$

$$
\textrm{for }-\infty < t <  \infty
$$

is periodic and its period $1/2$.

But I think it is not equal $ 1/2(\cos 4\pi t)  $ at $ t = 0 $ so it is not periodic because at $ t = 0 $ it is equal to $ 1 $ different from other multiples of period.

I have also curious about the value of $ u(t)+u(-t) $. Is it $0$, $1$, or $2$ ?

It is a matter of convention. The answer depends on the definition of the step function $u(t)$. In signal processing, it is common to define $u(0)=1/2$. Under this definition, we have that $Ev(\cos(4\pi t)u(t)) = 1/2(0.5+0.5)=1/2$, so the signal is indeed periodic.

If $u(0)$ is defined as either 0 or 1, then you'd be correct and the signal wouldn't be periodic.


Understanding Digital Audio - "From Audio Device" in Matlab I understand that digital audio has both a sample rate (time "discretization") and a bit depth ("amplitude" discretization). Also, it can have multiple channels (stereo would be 2 channels).

But things that I can't find readily in the internet:

What is the range that is discretized in the bit depth? My experience is in motions, so I know if I have an accelerometer with a range of +-10 g's, and I discretize it in 8bit, well I would have a resolution of 20g/256. What would be the equivalence in audio?
Matlab's Simulink has a block for reading audio. It's output has the dimensions of MxN(xTime), where M is the number of "consecutive samples" and N the number of channels. What is "consecutive samples" (also called in the same help file as "Frame Size (samples)"). If I set this value to one, I would get what I would expect from a Mono sound, just a time history of amplitudes. But if I set it to two, it's almost as if they are both the same with minor differences in fast transitions. Is this meant for to get a average of the signals to be somewhat more "true to sound"?

Here is the link for Simulink's block: http://www.mathworks.com/help/dsp/ref/fromaudiodevice.html

1) The microphone converts the audio signal to an electrical signal, and then that electrical signal often goes through an amplifier and then to the ADC.  The ADC usually samples the voltage of the electrical signal, which is directly analogous to the pressure of the audio signal.

If you want to think of it in audio terms, probably the most useful way to think of it is that your microphone/amp/something saturates at a certain dB level.  Let's call that "X" dB.  That is where your signal maxes out.  If you have 16-bit sampling then you can span a decibel range of $20log_{10}(2^{16})$ = 96 dB.  That means that, in an ideal world, you could cover audio signals from "X" dB to "X" - 96 dB.  Since the world is not ideal you would probably not get that entire range, but you would probably get pretty close.

2) "Consecutive samples" just means the number of samples.  They throw "consecutive" in there just to say that they are not skipping samples.  If you set that value to 1 you will only get one sample per channel, which is probably not what you want to do.  The number of channels refers to the number of audio sources.  For mono you would set that to 1.  For stereo you would set it to 2.


Compute the Second Order Derivative of Digital Image with Finite Differences I was looking for how to compute second order derivative of an image and came across the question kernels to Compute Second Order Derivative of Digital Image. In the top voted answer, it gives an example to get $3\times 3$ kernel $I_{xx}$ and $I_{xy}$. The 1st order derivative filter is defined as
$$
\left(
\begin{array}{cc}
 -1 & 1 \\
 -1 & 1 \\
\end{array}
\right)
$$
Then

$$
 I_{xx} = I_x \cdot I_x = \left(
\begin{array}{cc}
 -1 & 1 \\
 -1 & 1 \\
\end{array}
\right)\cdot\left(
\begin{array}{cc}
 -1 & 1 \\
 -1 & 1 \\
\end{array}
\right)=\left(
\begin{array}{ccc}
 1 & -2 & 1 \\
 2 & -4 & 2 \\
 1 & -2 & 1 \\
\end{array}
\right)
$$

and

$$
I_{xy} = I_y \cdot I_x = \left(
\begin{array}{cc}
 -1 & -1 \\
 1 & 1 \\
\end{array}
\right)\cdot\left(
\begin{array}{cc}
 -1 & 1 \\
 -1 & 1 \\
\end{array}
\right)=\left(
\begin{array}{ccc}
 1 & 0 & -1 \\
 0 & 0 & 0 \\
 -1 & 0 & 1 \\
\end{array}
\right)
$$

The computation of $I_{xx}$ makes sense to me, however, for $I_{xy}$, shouldn't it be defined as below?
$$
I_{xy} = I_x \cdot I_y = \left(
\begin{array}{cc}
 -1 & 1 \\
 -1 & 1 \\
\end{array}
\right)\cdot\left(
\begin{array}{cc}
 -1 & -1 \\
 1 & 1 \\
\end{array}
\right)=\left(
\begin{array}{ccc}
 -1 & 0 & 1 \\
 0 & 0 & 0 \\
 1 & 0 & -1 \\
\end{array}
\right)
$$

Pay attention that convolution mean flipping the kernel both on the x and y axis.
Hence the first element is a multiplication of $ -1 \cdot -1 $ which yields $ 1 $ as in the answer in the original post.
Pay attention that this is a discrete approximation of the gradient based on Finite Differences. This specific one is based on the Sobel 1st derivative filter.
In Finite Differences Coefficients you may find more variations.
As usual, the longer the filter the better the accuracy of the approximation with longer tarnsients.


Question about Fourier descriptor and the physical meaning of the coefficients I would have some questions about the Fourier descriptor and the physical interpretation of the individual coefficients. I'm rather new to this concept, so I would be looking for a simple answer.

To start with I have a curve in xy [mm] coordinates (188x2 points). To calculate the descriptor I use the function found in R.C. Gonzales 'Digital Image Processing using Matlab' (p. 629). From that I get a vector with the descriptors An (188x1).  As far as I understand the concept of the descriptor I can 'filter' out the high frequencies by reconstructing the shape again with only a handful of descriptors.

Here are my questions now:

I read that, the DC component A0 is the position of my shape in the image and its the mean values of the x,y coordinates. I understand that, I could calculate the position directly form my x,y coordinates, but I was wondering where this DC component is in the descriptor vector, since its a even number of values. Is it the term with the highest value?
Which brings me to my second question. If I recalculate the shape with only 2 descriptors I obtain a circle around the mean position with a certain radius. Is the radius then the second highest value?
Last question. Im mainly interested in the ellipses that is created with recalculating the shape with 4 factors. But I would also like to understand the whole concept correctly. I read that the A-1 and A1 coefficients represent the major (d) and minor (e) axis of this ellipses. And that with d = 2(|A-1| + |A1|) and e = 2(|A-1| - |A1|) the axis can be calculated in a length unit. Do I have to scale these values somehow? Because if I pick the third and fourth highest values, the result is not what I expect.

I hope someone can help me better understand this. Thank you in advance.

This question overlaps, to an extent, with this one. The best way to understand the Fourier Descriptors is through a basic understanding of the Discrete Fourier Transform.

The "curve" is basically a complex signal at the input of the DFT that is decomposed to a sum of wavy ellipsoid rings of increasing frequency.

Another convenient way to think about this is that what we usually see as "amplitude vs time" has now become "radius vs angle". Amplitude of a quantity versus time is how this quantity increases and decreases over time. In the radius vs angle way of thinking, you map the amplitude to the radius of an ellipsoid and the angle to a point on its circumference. Just as you read a diagram left to right, you "read" a curve beginning at some point and proceeding to the next either clockwise or counter clockwise.

  ...where this DC component is in the descriptor vector, since its a even number of values. Is it the term with the highest value?

The DC component is the first "bin" of the DFT. Given some signal $v[n]$, its DFT is $V = \mathcal{F}(v)$. But $V$ is usually indexed by $k$, as in $V[k]$. The $k^{th}$ value of $V$ is called a "bin". $V[0]$ is what is known as the DC component and it is basically a constant value. Because of the way the DFT is evaluated, the $V[0]$ is basically $\frac{1}{N}\sum_{v \in V} v$ and because $v$ is complex, the result can be complex too. But since you have encoded your $x,y$ curve as a complex signal of $v[n] = x[n] + j y[n]$, that $V[0]$ is basically representing the centroid of the curve.

  If I recalculate the shape with only 2 descriptors I obtain a circle around the mean position with a certain radius. Is the radius then the second highest value?

It is difficult to know what your reconstruction looks like without a specific example. The first component will be an ellipsoid which represents the mean radius through all points. If you keep adding components, all that you do is adding "squiggliness" to the average ellipsoid until a point where that original ellipsoid has become as squiggly as the curve.

Here is an excellent example:

From this paper.

Every value of $V$ is how quickly the radius varies versus the angle argument around the curve. So every value of $V$ is a "radius".

  I'm mainly interested in the ellipses that is created with recalculating the shape with 4 factors. But I would also like to understand the whole concept correctly. I read that the $A_{-1}$ and $A_1$ coefficients represent the major and minor axis of this ellipses. And that with $d = 2(|A_{-1}| + |A_{1}|) $ and $e = 2(|A_{-1}| - |A_1|)$ the axis can be calculated in a length unit. Do I have to scale these values somehow? Because if I pick the third and fourth highest values, the result is not what I expect.

What you are probably experiencing there is the scaling factor of the DFT.

I do not have Gonzalez' book handy right now but I will check and amend this part if it refers to something else, especially with respect to the way $A$ is indexed here.

In general there is a $\frac{1}{N}$ scaling factor for the DFT which is usually applied in the inverse transform (see for example here and here. So, the forward transform evaluates sums directly and this is why you end up with $V$s that are very high compared to the original signal. If you want to use those frequency coefficients directly, you have to divide them by $\frac{1}{N}$.

But, apart from the DC component, the rest of the ellipsoids do not have a "radius".

What are you trying to measure?

Hope this helps


Purely theoretical question about idelal filters and infinite oscillations I am asking perhaps a naive question, but still, it would be nice to have this formally stated one time:

in theory, if we could do it (of course, we cannot, but imagine we could), if we could convolve an upsampled image with an infinite sinc function (resulting from the ideal low-pass filter in freqency domain, i.e., the brick-wall) would the interpolation really be perfect, i.e., no Gibbs ripples ?

Note: historically the Gibbs phenomenon was described over a continuous-time domain periodic square wave which was constructed from a finite number of its Fourier series coefficients, however, one can interchange the domains and observe the same behaviour on the frequency response of a finite length discrete-time filter with its impulse response samples designating the Fourier series coefficients.

This question is essentially asking the consequences of the unavoidable Gibbs phenomenon on the mathematical construction of an ideal brick-wall filter whose frequency domain definition apparently (and confusingly) does not involve the expected Gibbs ripples...?

The practical answer is that as the filter length approaches infinity the peak magnitude of the ripples does not reduce, but its duration shrinks to a single point at the edges of the passband. From Riemannian integration point of view, this will not have an effect on the filter's output as if there's no Gibbs phenomenon.

On any finite length approximation of an ideal brickwall filter, however, the ripples do exist, and will alter the actual output from that of an ideal filter's computation. The better solution to reduce those ripples is to make use of a window instead of increasing the filter impulse response truncation length.


How to obtain the exact value of wavelength from a 2D FFT amplitude vs wavenumber plot like it is obtainable from 1D FFT amplitude vs wavenumber plot? I have a two dimensional multi modal spatial signal generated from a MATLAB code using sinusoidal functions of different wave numbers, amplitudes and phases. What I want to know is that if I have the amplitude vs wave number plot of that signal, how can I extract the wavelength of the different spatial structures. For 1D signal and 1D FFT I know that it is possible to extract the wavelength from the amplitude vs wavenumber plot by simply taking the reciprocals of the wave numbers with non zero amplitudes. But how can it be done for 2D signal and 2D FFT? I am attaching my code below.
Ws=160; % Sampling wavenumber @ 160 Hz
L=10; % Length of domain = 10cm
N = L*Ws; % Length of signal
x = (0:N-1)*(1/Ws); % Space vector
y = [(0:N-1)*(1/Ws)]'; % Space vector
x = repmat(x,1600,1); % Space matrix
y = repmat(y,1,1600); % Space matrix
V = sin((4pix)/L)+sin((4piy)/L); % Function in the spatial domain
`
fx = linspace((-Ws/2),Ws/2,N); % computing wavenumber vector fx
fy = [linspace((-Ws/2),Ws/2,N)]'; % computing wavenumber vector fy
fx = repmat(fx,1600,1); % computing wavenumber matrix fx
fy = repmat(fy,1,1600); % computing wavenumber matrix fy
T = fft2(V); % 2D FFT
subplot(1,2,1)
pcolor(x,y,V) %%%% Contour plot in spatial domain
ax = gca;
ax.LineWidth = 2;
colormap jet
colorbar
shading interp
pbaspect([1 1 1])
xlabel('X(cm)-->.')
ylabel('Y(cm)-->.')
title('V = sin((4pix)/L)+sin((4piy)/L)')
subplot(1,2,2)
pcolor(fx,fy,abs(fftshift(T))) % Contour plot of FFT amplitude vs wavenumber
lin = zeros(1,N);
hold on
plot(lin,fy(1:1600),'b') % Vertical line passing through (0,0) in fft amplitude plot
hold on
plot(fx(1:1,1:1600),lin,'k') % Horizontal line passing through (0,0) in fft amplitude plot
colormap jet
colorbar
shading interp
pbaspect([1 1 1])
xlabel('Wx(cm^{-1})')
ylabel('Wy(cm^{-1})')
title('FFT amplitude')
I am attaching the plot which I got as well. Any help would be greatly appreciated.

First, you want to find the peaks in this 2D plot. In your case, you know that there will be 2 peaks in your FFT (because you added two sines). Your plot shows 4 peaks because the FFT magnitude is symmetric, since your signal is real. You can discard negative wave numbers and look for peaks among the positive wave numbers.
The wave numbers corresponding to those peaks belong to your signal. There will be two wave numbers corresponding to each peak, $k_x$ and $k_y$. We know that $k = \frac{2\pi f}{c}$, where $k = \sqrt{k_x^2 + k_y^2}$ is the wave-number, $f$ is the frequency in Hz and $c$ is the speed of sound in the medium. Now, using $\lambda = \frac{c}{f}$, we can deduce that $\lambda = \frac{2 \pi}{k}$.


Definition of difference between compression and filtering I'd like to explain the motivation.

It sounds like a neuroscience topic but I am after the strict Signal Processing POV.

(I don't come from a Signal Processing background).

Interested to learn whether it's a rookie question in Signal Processing or  touches on a deeper paradox.

Motivation for Question
A scientist, many years ago, told me that the most advanced data compression we have is in the optic nerve. That visual data is compressed there and then bought to our brain for more processing.

(I'm guessing visual data is defined as a 3D field that maps light, based on a variety of factors e.g. intensity, wavelength, Cartesian coordinates etc).

If my dog and my cat and I  all look at my cup of coffee- presumably, we will all compress the visual data distinctly. For sure, we all use visual data very differently, which must mean we extract different meaning from it, ie. have different filters for what's relevant in the visual field (e.g. color, movement, sensitivity to small or large shapes).

Question

I am interested to learn whether, in Signal Processing, compression and filtration are separate concepts or they are strictly equivalent?

Because the more I thought about it, the less I can understand the difference between compression and filtering.

It seems like compression involves filtering automatically. (i.e. compression is unconscious filtering, whereas filtering is conscious).

Thanks!

My personal view: it depends a bit whom you talk to.

In the SP community, if people talk about "filters", more than 90% of the time, they will refer to linear filtering, in particular filters that can be described by $n$-th order linear differential equations, leading to rational functions in the Laplace/z-Domain. For such filters, we have an enormeous amount of tools to describe their input/output behavior and use this to design the filter coefficients with respect to all kinds of criteria.

On the other hand, when people talk about compression, it will in general involve nonlinear methods, since a purely linear compression is just not that efficient (sidenote: compressed sensing does exactly that, but it is not a compression method really, it's about the efficient sampling of analog signals and yet often mistaken for a compression method). After all, we would expect two very similar input signals to be compressed to the same representation and this already means there is some nonlinearity involved.

That said, some filters that are very commonly used are nonlinear, like the median filter, and there is no real reason to limit filters to linear ones. One could consider a much wider range of nonlinear filters. It's just in general harder to design these, at least from a classical filter design point of view.

You might still refer to a compression scheme as a filter and some people in SP will do that. I'd say it's a bit less common but that's just my point of view that is also heavily biased by my background.

Summarizing, also in response to your comments:

In compression, the goal is to represent a given signal with a minimum amount of data. If digital data is already available, the goal is thus typically to reduce the amount of data.
In filtering, the goal is typically to extract certain parts of the signal by suppressing others. This means that the amount of data may not change at all in a filtering process. Yet the term is so general that in some cases it might.
Compression will apply some sort of filtering (and maybe more).
The terms area not used consistently in the SP community.


wavelet transform boundary I am doing 1D wavelet decomposition and I am particularly interested in parts at the border of the signal. These parts are affected by boundary effects.
I know that some of the methods to extend the signal are periodic and reflexion, but none of these solutions are acceptable in my case.
I would like to extend the signal by the value at the border (something a bit like zero padding, but with a non zero value)

Does something like this has already been tested for wavelets?

To better extend the signal at the border, you can consider a constant extension: $[x_0,x_1,x_2]$ can be extended to the left (and similarly to the right) as $[x_0,x_0,x_0,x_1,x_2]$. However, this tends to break the slope, as it assumes a zero-derivative (flat signal).

A more advanced option resides in anti-symmetric extensions: they are sometimes called whole-sample or half-sample extensions, the two main options are:

$$[2x_0-x_2,2x_0-x_1,x_0,x_1,x_2]$$
and
$$[2x_0-x_2,2x_0-x_1,x_0,x_0,x_1,x_2]$$
In other words, the average of the extreme points is  $x_0$. This better respects derivatives in low noise conditions. Other types of extrapolation exist, more ad hoc (linear or parabolic fits), but may help. In the Matlab wavelet toolbox, Border Effects are described, and some Discrete wavelet transform extension modes are:

'asym' or 'asymh': Antisymmetric-padding (half-point): boundary value antisymmetric replication
'asymw': Antisymmetric-padding (whole-point): boundary value antisymmetric replication
'spd' or 'sp1': Smooth-padding of order 1 (first derivative interpolation at the edges)
'sp0': Smooth-padding of order 0 (constant extension at the edges)

As said by @DaBler, such extrapolations may require an increase in the transformed coefficients but:

for mere processing, you can crop the extended parts afterward,
there exists symmetric/antisymmetric filter-banks (e.g. $M$-band wavelets, $M >3$ and even) for which the extension can be made non-expansive.

Aside, if you really care about signals close to the border, you should take care of the sub-sampling, and possibly avoid it, as you cannot afford loosing samples. And more importantly, make sure that your wavelet is (close to) symmetric as well, unless you can lose the positive effect of your extension.

Finally, wavelets build on the interval adapt the wavelet size close to the border.


How to transmit the real+imaginary output of an IFFT over a transmission channel? I'm working on a modem system using a 16-QAM scheme. The transmitter runs an IFFT to modulate the signal, and has a 12-bit DAC. The receiver runs an FFT to demodulate the signal, and has a 12-bit ADC.

Here are the general parameters:

Sample rate is 512,000 SPS
FFT size is 64 samples
4 bits per bin/frequency with a total of 10 bins/frequencies.
Frequencies are in the range 32KHz to 104KHz. Frequency spacing is 8KHz.
The rectangular 16-QAM constellation diagram is similar to this: http://en.wikipedia.org/wiki/File:QAM16_Demonstration.gif

The IFFT on the transmitter produces two 64 element arrays, real and imaginary. The transmitter has a 12-bit DAC available to transmit the signal.

I'm having trouble understanding exactly how to transmit this output of the IFFT? I've read in theory that I need to combine the two arrays (imaginary and real parts) via sine and cosine, but am having trouble understanding how this concept practically gets done in terms of programming.

Can anyone outline briefly (maybe a snippet of code or something?) how I combine the two 64 element arrays of real and imaginary data to feed into the DAC to transmit the signal on the wire?

If you want to transmit a complex signal, you need two independent channels: one for the real part of the signal and another for the imaginary part. This makes sense, as the complex plane is two-dimensional. The question is, how can those channels be realized? It depends on what equipment you have available.

If you have multiple wires/DACs that can be used to convey the signal, then you can just use baseband transmission. Route the real part of your signal to one DAC and the imaginary part to the other. These would then feed separate transmission lines to carry the signals to the receiver. This method isn't often used, though. Why use multiple wires when you can just...

If you have a single wire/DAC that can be used to transmit the signal, then you can use carrier modulation. Basically, you would modulate the real and imaginary parts of the signal onto separate sinusoidal carriers (this is done by multiplication by the desired carrier), which are then observed at the receiver and reassembled into the composite complex signal that you started with.
A common technique for modulating complex signals is called I/Q modulation. This scheme uses a pair of sinusoids at the same frequency that have a phase offset of 90 degrees. This particular phase offset results in carrier signals that are orthogonal to one another, which allows the receiver to (in a perfect world) losslessly recover the two modulated signals. In practice, there will always be some "crosstalk" between the I and Q channels (which you could think of as carrying the real and imaginary parts of your modulating signal, respectively).

In your question, you suggested that you have a single DAC, so I would go the carrier modulation route. You'll need to select a carrier frequency appropriate for your application, taking the following things into account:

The bandwidth of your modulating signal.

The propagation medium to the receiver. If this is a wire, then you can often get by with as low of a carrier frequency as your signal bandwidth will allow (barring any limitations imposed by the path to the receiver, such as transformer coupling or other details). Wireless links often have quite high carrier frequencies (in the tens, hundreds, or thousands of MHz), depending upon characteristics such as desired range, transmit power, licensing, propagation environment, etc.

Hardware limitations in your transmitter/receiver. You need to generate the carrier-modulated signal somehow! If you don't have analog hardware to do that for you (i.e. you output a pair of baseband signals and an analog transmitter upconverts it to an I/Q carrier-modulated output), then you'll need to synthesize the carrier-modulated signal digitally. That may or may not be possible given your system's constraints on DAC sample rate, processing throughput, and so on.

FYI, the design that you described for your system sounds a lot like orthogonal frequency-division multiplexing (OFDM), albeit with a lot of unused carriers. You might find some useful information if you study up on how those systems are usually implemented.


Power of periodic signal Given a signal
$s(t) = e^{j(t+\pi)}$,
I conclude that the signal is periodic with period $T=2\pi$, so its power should be
$$P = \frac{1}{T}\int_{-T/2}^{T/2}|s(t)|^2dt= \frac{1}{2\pi}\int_{-\pi}^{\pi}|1|^2dt = \frac{1}{2\pi}t|^{\pi}_{-\pi}=\frac{\pi-(-\pi)}{2\pi}=1$$

Is the formula for the power of a periodic signal correct?

The average power of a signal $s(t)$ is given by

$$\overline{s^2(t)}=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{-T/2}^{T/2}|s(t)|^2dt\tag{1}$$

If $s(t)$ is periodic, $(1)$ is equivalent to the average power in one period:

$$\overline{s^2(t)}=\frac{1}{T}\int_{-T/2}^{T/2}|s(t)|^2dt\tag{2}$$

where $T$ now denotes the period of $s(t)$.


Is ripples in frequency response of a filter ,result of periodic extension of a signal in time domain thereby producing discontinuities? I am aware that due to gibbs phenomenon, ripples are induced at the vicinity of transition in the frequency response of a filter and ripples occur in frequency domain when there are discontinuities in the time domain signal.

Now my doubt is how and why discontinuities are present in time domain signal?

Is it due to periodic extension of aperiodic signal in time domain to compute DFT?

please correct me if I am wrong in my understanding and provide an explanation.

Thanks!!

It's mostly related to what's called the uncertainty principle of the Fourier Transform: "the more compact a signal is in one domain the more spread out it is in the either". A sharp transients in the time domain has a wide spectrum. A sharp transition in the frequency domain has a long time domain ringing. See for example http://www.ams.org/samplings/feature-column/fcarc-uncertainty


How to calculate step and impulse response at multiple sampling rates I'm struggling with a basic concept.  Suppose I have a system $y(t)=e^{-at}u(t)\ast x(t)$ and I want to plot impulse response and step response using Matlab conv() function at various sampling rates.  I'm having trouble understanding what gain I should apply and why.  I tried:

a = 0.05;
cutoffTime = 250;
T = 1/samplingRates(i);  %assume this is defined somewhere [1, 10, 20]
timeRange = 0:T:cutoffTime;
step = ones(1, length(timeRange));
impulse = zeros(1, length(timeRange));
impulse(1) = 1;
exponential = exp(-a*timeRange);
stepResponse = conv(exponential, step, 'full');
impulseResponse = conv(exponential, step, 'full');

I understand that mathematically sampling is modeled as multiplication by a pulse train, which introduces gain of $1/T$ to the discrete system.  So the reconstruction filter needs gain of T.  In the above code, for different sampling rates (say 1, 10, 20 samples/second) if I don't apply the gain of T to step response, I get different sized step response.  However if I apply the gain T to impulse response, I don't get the impulse response I expect for all 3 sampling rates.  Can anyone clarify where I'm going wrong in my implementation?

Thank you.

This is not the full answer you want, but the problem is that you can't permute sampling and convolution. Let's say that the $D(\cdot)$ operator takes a continuous signal and returns its sampled version. What you're doing is

$$y'(n) = D\left(e^{-at}u(t)\right) * D\left(x(t)\right)$$

What you want to get is

$$D\left(y(t)\right) = D\left(e^{-at}u(t) * x(t)\right) \neq y'(n).$$

It is in general impossible to do what you're trying to do, but I'm sure that good approximations exist. I'm not an expert on what they are or how to use them.


Deblurring algorithm to precede thresholding - speed over accuracy I'm writing an app that recognizes Sudoku puzzles from a camera input.  I'd like to remove camera blur from the images to improve recognition.  Here is an example image:

Since I'm processing a continuous camera feed, I'm concerned about speed, since slow algorithms will degrade my frame rate.  On the other hand, the technique doesn't have to be super-accurate, since I will be thresholding/binarizing the image afterwards.  The deblurring algorithm I use only needs to be accurate enough so that the digits will be clear after thresholding.  Here is what the image currently looks like after a segmented Otsu threshold (the threshold level is selected independently for each 3x3 grid of cells):

I'm looking to write an implementation that's compatible with OpenCV data structures.  This means that an off-the-shelf Matlab solution isn't an option.  However, I found the Matlab documentation on image deblurring to be quite helpful in understanding the general principles, as well as this site.  The Matlab functions are all computationally-intensive in any case, and I wasn't able to pick an estimated point-spread function (PSR) that gave good results.

What are your suggestions for a fast deblurring algorithm that allows features to be distinguished easily by a subsequent threshold?

You only need to find and report positions of nine specific symbols.  You can do this by computing correlations (done quickly by multiplying in Fourier space) with samples of the nine digits snipped from a typical puzzle image, blurred as it is.  All blurred '9' will look the same, and likewise all other digits.  I'm assuming puzzle images are taken with the same camera the same way.

It might work better working on high-pass filtered versions of the input puzzle image and image snips, by subtracting from each a more-blurred copy.


Phase shift in higher order filters I was doing some research on filters, So I designed butterworth filter based on pole locations. I applied to a simple sine signal to check for the outputs. I cannot find the reason for what I noticed.

Why does higher order filter have alot of phase shifts. Is it trying to adjust for more noise reduction? If, Yes how is it mathematically possible? I have only changed the order of the filter, I have kept cut off frequency same

Discrete filters are divided in two major classes: Finite Impulse Response (FIR) filters and Infinite Impulse Response Filters (IIR).

Both types of filters are structured around a fundamental unit. A fundamental operation that is responsible for the filtering. That is, the delay unit. Of course, besides the delay unit you also have the sumation points. But sumation, on its own, without delay doesn't produce filtering.

So, both filters, produce outputs that depend on the past of the inputs (and / or outputs).

The concept of "filter order" is defined for both types of filters but in a slightly different way (it's still the same thing really).

Let's look at an FIR filter:

$$y_n = \sum_{m=0}^{Nh-1} x_{n-m} \cdot h_m$$

Here, to produce one value for $y_n$ we sum $m$ products of the past of the signal and some impulse response $h$. $Nh$ is the total length of $h$.

The order of the filter here is $N_h$, directly proportional to the number of delay units in the filter.

So, by increasing the order, we do get a sharper filter but also a longer filter and the longer the filter is, the longer it takes for one of the $x_n$ to propagate through the series of delays.

And that is where the increased group delay is coming from! Furthermore, in an FIR filter all frequencies are delayed proportionally. So, bass frequencies are delayed less than treble frequencies in a predictable (and linear) way because the time it takes for the sinusoids to propagate through the delay line is fixed.

Now, let's look at an IIR filter:

$$y_n = \sum_{m=0}^{Nh-1} x_{n-m} \cdot h_m - \sum_{k=1}^{Ng-1}y_{n-k} \cdot g_k$$

Which, if you take a closer look, is like two FIR filters back-to-back. The left part, we have seen before and it is the FIR filter from the first equation. The right part is an FIR but this time it is on past values of the output. We therefore have introduced feedback...and all the "joys" that come with it.

The order of an IIR filter is now the length of the feedback delay line (i.e. the $Ng$ here) and because of this interaction between the input and output, the group delay of an IIR filter varies with frequency. So, treble frequencies might be going through much faster than bass frequencies in an IIR filter.

OK, so now let us look at your experiment. You design a filter and you pass through it a simple sinusoid at some fixed frequency. Then you increase the order and you observe a longer delay.

Which is perfectly reasonable, because, irrespectively of the type of filter (FIR / IIR), increasing the order, means adding more delay units.

Hope this helps.


Bit scrambling as last state of transmitter - what advantages are there? I believe that I read, for digital modulation with multiple levels, that the AGC makes use of the fact that the symbol sequence, before noise etc., can be modeled as a uniform distribution.

$$\text{Source $\rightarrow$source-coding $\rightarrow$FEC$\rightarrow$ interleaving $\rightarrow$ bit-scrambling $\rightarrow$ modulation $\rightarrow$ channel}$$

What other advantages are there to bit scrambling (not for privacy/security).

Chris.

In normal usage, scrambling is used for two reasons:

It is used to remove the possibility of a long sequence of 1's and 0's in the bit sequence. The long sequence of 1's and 0's make timing synchronisation and clock synchronisation tougher at the receiver as regular transitions  help in working  of adaptive circuits like AGC and phase locked loop.
It eliminates the dependence of signal's power spectrum on the transmitted information sequence thereby keeping it below the maximum power spectral density requirement. If scrambling is not done, power might be concentrated in a narrow frequency band thereby causing intermodulation and crossmodulation distortion to adjacent channels.

Reference:
http://www.princeton.edu/~achaney/tmve/wiki100k/docs/Scrambler.html


Convert Fourier transform of a signal to $\omega$ form? Suppose we have this signal:

$$x(t)=5 \sin(2\pi 1000 t) \cos(2\pi 10000 t)$$

I calculated the Fourier transform as below:

$$X(f)=\left(j\frac{5}{4}\right)\left[\delta(f-9000)+\delta(f+11000)-\delta(f-11000)-\delta(f+9000)\right]$$

and then I would like to convert it to $\omega$ form. I know that the result is

$$X(\omega )=\left(j\frac{5}{2}\right)\left[\delta(\omega -9000)+\delta(\omega +11000)-\delta(\omega -11000)-\delta(\omega +9000)\right]$$

I know that
$\omega=2\pi f T$
where $T$ is sampling rate.

What I cannot understand is how we convert the Fourier transform of our signal to $\omega$ form using above formula.

In your example there is no sampling involved, so you simply have $\omega=2\pi f$, where $f$ is the frequency in Hertz, and $\omega$ is the frequency in radians.

Your first result for $X(f)$ looks correct. However, your result for $X(\omega)$ is wrong. What you need to know to be able to convert an expression with Dirac deltas from $f$ to $\omega$ is the following relation:

$$\delta(af)=\frac{1}{|a|}\delta(f)\tag{1}$$

With $a=2\pi$ and $\omega=2\pi f$ you get from $(1)$

$$2\pi\delta(\omega)=\delta(f)\tag{2}$$

You get the two following basic Fourier transform correspondences, depending on whether you use $f$ or $\omega=2\pi f$ as the frequency domain variable:

$$\begin{align}e^{j2\pi f_0t}&\longleftrightarrow \delta(f-f_0)\\
e^{j\omega_0t}&\longleftrightarrow 2\pi\delta(\omega-\omega_0)\end{align}$$

So your expression for $X(\omega)$ misses a factor of $\pi$, and the sums and differences of the two frequencies in the arguments of the Dirac deltas need to be multiplied by $2\pi$:

$$X(\omega)=\frac{5\pi j}{2}[\delta(\omega-2\pi9000)-\delta(\omega+2\pi9000)+\delta(\omega+2\pi11000)-\delta(\omega-2\pi11000)]$$


Multisampling of an ADC channel in a PWM cycle in TI C2000 I want to sample an ADC channel multiple times in one PWM cycle.
I can trig any ADC channel in anywhere in my PWM timer counter now. But I want to trig and sample same channel in many times within a pwm period.
I tried to explain by drawing below.

Is there any solution for this?
Thanks and best regards

Assuming your PWM above is ePWM1.
1 - Configure ePWM2 to be synchronized to ePWM1.
2 - Adjust ePWM2 period to be N times smaller.
3 - Now trigger the ADC sampling to ePWM2. Assuming you sample on the bottom and top of ePWM2 triangle, you will sample twice per ePWM2 period and 2*N times per ePWM1 period.


EEG Sampling, Which Order should cover my result? I have EEG data sampling rate is 250 Hz and I want to band pass this data in ( 0.5-3 Hz &  3-8 Hz &  8-13 Hz & 13-30 Hz & 30-40 Hz).

Currently, the data is in the frequency range of 0.5-40 Hz. (my implementation code is in Matlab, filter designer and EEGlab).

Is the Butterworth filter good here?
How do I apply Butterworth filter on this data?
What should be the ideal filter order for get the better result?

Well, you're looking for a filter that can produce steep flanks, and suppresses power well over all the stopband. Sounds like a job for a DFT approx/window filter design with a Dolph-Chebychev window.

In matlab, you'd use the fir1 function to design your bandpasses with the chebwin giving you an appropriate windowing function.

You'll need to specify an order as first argument to fir1.

Since we're not really computationally constrained, let's use the same order for all band passes, as that will lead to constant group delay.

What you forgot to specify is how much overlap you can tolerate: no real filter in this world can cut off at e.g. 3 Hz sharply! You need to allow for some transition width. That transition width defines the order you need (the sharper, the higher the order.

Following this answer, assuming 0.1 passband ripple and 60 dB stopband attenuation, we'd get for a transition width of 0.1 Hz for the first filter an order of about 1000.

That's a large filter, and it will probably take a while for the fir1 method to complete, but because you have so little data to work with (in terms of computer speeds), that's going to be alright when applying the filter using the filter function.


2D Fourier Transform of Rotated Discrete Domain Signal Assume we know that the Fourier transform of a signal $x(n_1,n_2)$ is $\mathcal{F}(x(n_1,n_2))=X(\omega_1,\omega_2)$. What is the Fourier transform of the signal after being transformed by a rotation matrix?

I have found a property of the DFT that states that rotating the spatial domain contents rotates the frequency domain contents which is called the Rotation Property.  I am assuming this extends to DDFT's (Discrete Domain Fourier Transforms).  However, I would like to derive it.

We can represent our signal using vector notation as $x(\boldsymbol{n})$, where $\boldsymbol{n} = \begin{bmatrix} n_1 \\n_2\end{bmatrix}$.

The rotation of $\boldsymbol{n}$, denoted $\boldsymbol{n}' =\begin{bmatrix} n_1' \\n_2'\end{bmatrix}$, can be expressed as
$$\boldsymbol{n}' =R\boldsymbol{n} = \begin{bmatrix}
       \cos\theta & -\sin\theta \\[0.3em]
       \sin\theta & \cos\theta
     \end{bmatrix} \begin{bmatrix} n_1 \\n_2\end{bmatrix}
$$

So I am trying to find the Fourier transform of $x(R\boldsymbol{n}$).

My attempt is as follows:

Step 1: Write the Fourier transform in vector notation.

Define $\boldsymbol{\omega} = \begin{bmatrix} \omega_1 \\ \omega_2\end{bmatrix}$.  We can now write the 2D Discrete Domain Fourier transform as: $$\mathcal{F}(x(\boldsymbol{n}))=X(\omega_1,\omega_2)=X(\boldsymbol{\omega})=\sum_\boldsymbol{n}x(\boldsymbol{n})e^{j\boldsymbol{\omega}^{T}\boldsymbol{n}}$$ because $\boldsymbol{\omega}^{T}\boldsymbol{n}=\omega_1 n_1 + \omega_2 n_2$.
Step 2: Solve for $\boldsymbol{n}$:
\begin{align}
\boldsymbol{n}' &=R\boldsymbol{n}\\
R^{-1}\boldsymbol{n}' &=R^{-1}R\boldsymbol{n}\\
R^{-1}\boldsymbol{n}' &=\boldsymbol{n}\\
R^{T}\boldsymbol{n}' &=\boldsymbol{n}
\end{align}
Noting that $R^{-1}=R^{T}$ for rotation matrices.
Step 3: Apply Fourier transform to rotated signal:
$$\mathcal{F}(x(\boldsymbol{n}'))=\sum x(\boldsymbol{n'})e^{j\boldsymbol{\omega}^{T}\boldsymbol{n}}$$
Plug in $\boldsymbol{n} = R^{T}\boldsymbol{n}'$
$$\mathcal{F}(x(\boldsymbol{n}'))=\sum x(\boldsymbol{n'})e^{j\boldsymbol{\omega}^{T}R^{T}\boldsymbol{n}'}=\sum x(\boldsymbol{n'})e^{j(R\boldsymbol{\omega})^{T}\boldsymbol{n}'}$$
And therefore
$$\mathcal{F}(x(R\boldsymbol{n}))=X(R\boldsymbol{\omega})$$

Which I am interpreting as meaning the original statement of rotating the spatial domain contents rotates the frequency domain contents.

Can anyone confirm this is correct?

The Rotation Property means that the choice of coordinate direction would not affect the spectrum of the signal itself.

To perform the rotation, I think it is better to do in the spatial domain since frequency domain is a complex one and it will just increase the complexity.