Engine for Google / Browsers to detect manipulated images? Would it be possible for Google and web browsers to scan images and detect doctored photos?

For example, a web browser could check for abnormal color gradients and signals that indicate if a photo has been retouched?

Google could fingerprint photos into ASCII code, then compare to elements copied from older photos. If that was done in conjuction with people appearing in news articles, accidents, or porn, forgetting all the animal and technology photos, it wouldn't be a programming impossibility?

I was imagining a google images page where all the fake photos were shown.

Overview:
Detecting fake images is part of the research field of media forensics. This is an ongoing research niche.
There already exists a fairly mature set of computational tools that can identify various types of image manipulations. They can also indicate the region in the image that has been altered. Concurrently, the improvements in photo-editing tools, like Adobe Photoshop, had made the detection of manipulated images more challenging.

Would it be possible for Google and web browsers to scan images and photos, to detect fakes?

Yes, it is indeed possible. Integration into a popular web browser is probably pending, but the technology already exists. You can try it for yourself here
There is also an open source repository for those interested to contributed or modify here

For example, a web browser could check for abnormal color gradients and signals that indicate if a photo is fake?

A popular example of fake photo is the missile launch by Iran:

The actual lunch had only 3 missiles and 1 missile-launcher failed:

The manipulated image used copy-move forgery from the original image to synthesize the fake 4th missile. A clone-tool detection algorithm can detect this:

Google could fingerprint photos into ASCII code, then compare to elements copied from older photos. If that was done in conjuction with people appearing in news articles, accidents, or porn, forgetting all the animal and technology photos, it wouldn't be a programming impossibility?

There are examples of easy to detect fakes when the manipulator has used existing images. Typically this happens when the manipulator uses stock photography. So the answer is Yes, you are right! However, this is a small subset of fake images and more skilled manipulator have better ways of covering their tracks. Additionally, it does not work when the fake image is completely new.

I was imagining a google images page where all the fake photos were shown.

This question raises the point of provenance of images. The image in raw format from the sensor (camera) is the original. Every operation subsequently done to the image is by definition a manipulation. This includes saving the image in lossy format like jpeg. (Skilled image manipulators will cover their tracks by doing a double jpeg compression!) However, some manipulations are malicious, whereas others are not. There can indeed be a Google page with all fake photos (ie. photos with lower provenance than the 'original' photo). This might even become a requirement in photo-journalism in the near future!
There were multiple interesting sub-questions raised here. I hope I have at least partially answered most of them.


Fitting fixed number of complex exponentials to complex signal I have a complex data series, and I would like to fit a fixed number (in this case, two) of complex exponentials of the form $Ae^{jBn}$, where A, B are complex. Not interested in the phase (i.e. arg of A), just its magnitude, and mainly the complex frequency. Target length is maybe 100 samples.

I would like to be able to capture the periodic destructive interference as the addition of these two fitted complex exponential signals.

Complex time-series in 3D:

I have looked into Prony-type methods and noise subspace methods, but I'm at a loss as to what is applicable here.

Complex A and B are a bit messy. What about using the signal model

$$ f(t) = \sum_{k=1}^{2} A_k\, e^{j(w_k t + \phi_k)} $$

where $w$ is the frequency and $\phi$ is phase.

Since sinusoids are orthogonal, probably the easiest thing to do is to Fourier transform the signal and pick out the two largest peaks in the positive part of the spectrum. These will correspond to the two main sinusoids.


Proof of properties of Fourier series in CT I feel problem in understanding the proof of Fourier series properties

Time scaling
\begin{align}
b_k &= \frac{1}{T}\int_{T}x(t)e^{jk\omega_0t}dt\\
& = \frac{a}{T}\int_{T/a}x(at)e^{jk(\omega_0a)t}dt\tag{$\scriptstyle{a - \text{scaling factor}}$}\\
& = \frac{a}{T}\int_{T}x(L)e^{jk\omega_0L}\frac{dL}{a}\tag{$\scriptstyle{L=at}$}\\
&= \frac{1}{T}\int_{T}x(L)e^{jk\omega_0L}dL = a_k\\
\end{align}
For this I want to ask that how $T/a$ term in the integral changes to $T$, and $dt$ changes to $dL/a$?
Multiplication
\begin{align}
x(t)y(t)&= \sum^{+\infty}_{k\ =\ -\infty} a_k  e^{jk\omega_0t}\sum^{+\infty}_{k\ =\ -\infty} b_k  e^{jl\omega_0t}\\ &=\sum^{+\infty}_{k\ =\ -\infty} \sum^{+\infty}_{l\ =\ -\infty} a_k b_le^{j(k+l)\omega_0t}\\
&=\sum^{+\infty}_{m\ =\ -\infty} \left[\sum^{+\infty}_{l\ =\ -\infty} a_{m-l} b_l\right]e^{jm\omega_0t}\tag{$\scriptstyle{m=k+l}$}
\end{align}
I couldn't understand this last step, why the summation changes from $k$ to only  $m$? Shouldn't it be $m-l$ and also is it done here to prove the multiplication property? or there are more steps?

Looking at the change of terms in the equality:
$$\frac{a}{T}\int_{T/a}x(at)e^{jk(\omega_0a)t}dt = \frac{a}{T}\int_{T}x(L)e^{jk\omega_0L}\frac{dL}{a}$$
$$at = L\implies \begin{cases}\text{if }t=T/a\implies L=T\\adt=dL\implies dt=\frac{dL}{a}\end{cases}$$
Looking at the equality again:
$$ \sum^{+\infty}_{k\ =\ -\infty} \sum^{+\infty}_{l\ =\ -\infty} a_k b_le^{j(k+l)\omega_0t}=\sum^{+\infty}_{m\ =\ -\infty} \left[\sum^{+\infty}_{l\ =\ -\infty} a_{m-l} b_l\right]e^{jm\omega_0t} $$
$$m=k+l\implies \begin{cases}k=m-l\quad({\scriptstyle{\rm indeed}})\\\text{and } m-l=\pm \infty \implies  m=\pm\infty+l=\pm\infty\end{cases}$$


Windowing function squared? Reading this paper about WSOLA (Verhelst, Roelands) (PDF)

At page 2 there's a "time-warped STFT's synthesis function" involving $v(n)=w^2(n)$ which is a windowing function. A windowing function squared.

What's the interpretation/meaning of windowing function squared? What does it mean visually?

it's just another window.  usually this happens when a window is applied on the time-domain data going into the STFT (which is a DFT), some process happens to the frequency-domain data after the STFT, then there is an inverse DFT bringing it back to the time domain, and the window is applied again to that output time-domain data.  if that process done to the data in the frequency domain happens to be the null process, then what comes out of the iDFT is exactly what went into the DFT.  so then, in that case, the window is applied twice and that is the same as applying the square of the window once.

i have been having a little argument with someone on comp.dsp about this.  there is a term i hadn't seen before ("Princen-Bradley condition") which i think means the same as a "complementary-power" window.  if the square of a complementary-power window is applied to an STFT system above with a null process, the output will be exactly the input when 50% overlap is used.  but if the square of a complementary window (such as the Hann) is used, then, with 50% overlap you do not get perfect reconstruction.  so the square-root of a Hann window would be a good candidate for a complementary-power window.  and you would have to apply that window twice to the time-domain data to get perfect reconstruction.


What is exactly Power Spectrum Density? In papers or textbooks, I have found several different definitions of PSD. I think I understand the meaning of PSD and all these papers or textbooks agrees that PSD and Fourier Series share a proportional relationship. But still there are different formulas from P(f)=1/T(〖F(f)〗^2), P(f)=〖∆t〗^2/T(〖F(f)〗^2), P(f)=1/4∆f(〖F(f)〗^2), P(f)=(〖F(f)〗^2) or others. From what I learn, I think on the point of conception of power(or energy) and Parseval's theorem, discrete PSD should be described as simply P(f)=1/T(〖F(f)〗^2). So why is there so many different formulas? Are they wrong or do I misunderstand something?

Power Spectrum Density (PSD) represents the spectral distribution of total power of a given WSS random signal. Conventionally, it's defined as the Fourier transform of the Auto-Correlation function (sequence) of the WSS random process as:

$$ P_x(\omega) = \mathcal{ F } \{ r_{xx}[m] \} $$

where the auto-correlation sequence (ACS) $r_{xx}[m]$ of the WSS $x[n]$ is defined by:

$$ r_{xx}[m] = E\{ x[n]x^*[n-m] \} $$

These are theoretical entities, you may not compute them in practice where only a single or a few realizations of the random process are available. In such a case, the PSD is estimated and there are different estimation techniques, one of the most common and simplest is known as the periodogram estimate of PSD and given by

$$ \hat{P}_{x}(\omega) = \frac{1}{N} |X(\omega)|^2 $$

where $X(\omega)$ is the discrete-time Fourier transform of the observed signal $x[n]$ of length $N$.


should this be viewed as an impulse response or step response I'm trying to teach myself the relation between simple discrete ODE's and the impulse response-step response concept.

Getting back to the question: I don't expect anyone to read the whole thing but I'm confused early on which is a good thing.

Right at the beginning of the link, the author provides an example of a leaky integrator model in which the water facet is turned on and the water goes into the bucket and there is no hole in the bucket.

Given the leaky integrator model, the author then shows how, if there is
no hole in the bucket, then by setting k = infinity in the model, one can model the water level for the bucket. He then provides

a plot of the input: the water going into the bucket and
a plot of the output: the water level of the bucket.

My confusion stems from the plot of the input. I think that the plot of the input, rather than being a spike as is shown at time t = 1, should instead be a step response where there are spikes all the way from zero to one ? so basically a step response ?

My thinking is that, if the water is turned out at time zero, then the spike for the input  is at t = 0 and remains a spike until time t = 1 ?

The output plot seems fine to me.

If I am correct about the input plot being incorrect as is, then it kind of tells me that I understand the basic impulse response-step response concept. If I'm not, then, even though I loved the example, I'm still not quite getting it.

I asked the creator of the link in an email and he-she hasn't gotten back to me yet and I would like to understand this and don't honestly know if I'll get a response ( no pun intended ).

Rather than spending an efford on the link-exploration, I would instead state here the simplest and the most basic description of the impulse response and step response of causal LTI systems within the particularity of a continuous time leaky integrator. I assume that you know about linear time invariant (LTI) systems and basic differential equations.

A (continuous-time) LTI system is characterised by its impulse response $h(t)$ defined via the relationship $$h(t) = \mathcal{T} \{ \delta(t) \}$$ such that
$$ \delta(t) \rightarrow \boxed{ h(t) } \rightarrow h(t)$$

Where $\delta(t)$ is the impulse that is at the input. The step response can be defined as $$s(t) = \mathcal{T} \{ u(t) \}$$ such that
$$ u(t) \rightarrow \boxed{ h(t) } \rightarrow s(t)$$ where $u(t)$ is the unit-step function.

One can further relate the step response to impulse response as $$ s(t) = \int_{-\infty}^{t} h(\tau) d\tau \longleftrightarrow h(t)=s'(t)$$

Finally, for an arbitrary input $x(t)$ the corresponding output $y(t)$ of any LTI system is given by the convolution integral as:
$$ y(t) = h(t) \star x(t) = \int_{-\infty}^{\infty} h(t-\tau) x(\tau) d\tau = \int_{-\infty}^{\infty} x(t-\tau) h(\tau) d\tau = x(t) \star h(t) $$

Now before coming to the leaky-integrator, we shall consider the ideal integrator:

The ideal integrar is the system with the input-output relationship as
$$ y(t) = \mathcal{T}\{ x(t) \} = \int_{-\infty}^{t} x(\tau) d\tau $$

Considering that ideal-integrator is an LTI system we can regard the above integral as a convolution operator and hence deduce that the impulse response of the ideal integrator is $$h(t) = u(t)$$ then the step response will be $$s(t) = \int_{-\infty}^t u(\tau) d\tau = t u(t) = r(t)$$ $r(t)$ being the ramp function.

Note that the differential equation associated with the ideal-integrator can be obtained by diferentiating its input-output relationship with respect to time as

$$ \frac{d y(t) }{dt} = \frac{d}{dt} \int_{-\infty}^t x(\tau)d\tau $$ which results in (after following the differentiation of the integral)
$$ y'(t) = x(t) $$ with properly chosen initial conditions (according to the applied input) so that the ODE (an LCCDE) represents an LTI system.

The leaky-integrator can be defined in a number of ways as the term leaky may accept multiple mathematical definitions but one of the most common has the following input-output relationship:

$$y(t) = \mathcal{T} \{x(t) \} = \int_{-\infty}^{t} e^{-\alpha (t-\tau)} x(\tau) d\tau $$ Which has the interpretation that, when $\alpha > 0$, past values of the input are ignored more and more as time goes on, when computing the current value of the output. This is achived by the decaying-exponential weighting applied to the integrator input $x(t)$.

Again considering that this is an LTI system and therefore that its output is given by a convolution integral, then we can deduce that the impulse response $h_l(t)$ of the leaky integrator is
$$h_l(t) = e^{-\alpha t} u(t)$$ and the associated step response can be obtained as $$s_l(t) = \int_{-\infty}^t h_l(\tau)d\tau = \int_{-\infty}^t e^{-\alpha \tau} u(\tau) d\tau = \int_{0}^t e^{-\alpha \tau} d\tau = \frac{1}{\alpha} (1-e^{-\alpha t})$$

The associated ordinary differential equation of the leaky-integrator can be obtained similarly by diferentiating its input-output relationship with respect to time $t$ as

$$ \frac{d y(t) }{dt} = \frac{d}{dt} \int_{-\infty}^{t} e^{-\alpha (t-\tau)} x(\tau) d\tau .$$ After following the differentiation of the integral carefully one can reach
$$ y'(t) = e^{-\alpha (t-t)}x(t) t' - e^{-\alpha(t+\infty)}x(-\infty) \frac{d}{dt}(-\infty) + \int_{-\infty}^{t} -\alpha e^{-\alpha (t-\tau)} x(\tau) d\tau  $$
It can be seen that the first term is $x(t)$, the second term is zero and the last term is $-\alpha y(t)$, hence  we deduce
$$y'(t) = x(t) - \alpha y(t) \longleftrightarrow y'(t)+\alpha y(t) = x(t)$$ as the ODE that define the LTI leaky integrator.


Carrier frequency offset estimation for QPSK and asymmetrical spectrum I am working on estimation of carrier frequency offset for QPSK data from a satellite. From my understanding, there are two main approaches: either remove data from the carrier (by taking the fourth power), or with band-edge filtering, or some other method of measuring difference of energy content in positive and negative parts of the spectrum.

In my particular case I have the signal converted to base-band (and with some offset). I'm now trying to estimate that offset. I can correctly estimate offset by taking the fourth power of the signal, and then running an FFT on that signal. However, I'm struggling with band-edge filtering, and it seems that the issue is not my implementation of the filter, but rather, the signal itself. While the algorithm works like charm on synthetic QPSK data, when I run it on some real satellite data, it is not able to compensate for frequency offset. I tried to find the cause of this behavior, and I was able to find that the signal spectrum is asymmetrical, and negative frequency content is having more energy than positive, as shown in the figure below (figure shows FFT of signal with frequency offset removed). What I do not understand why there is such a difference in the energy in the positive and negative parts of the spectrum, and I would appreciate if somebody could provide me some insight into this.

A likely source of the amplitude imbalance given it is on the actual signal is filtering distortion along the signal chain such that the gain is not flat across the bandwidth of the spectrum. This would not be uncommon with analog filtering and can be compensated for with equalization (if needed).

For an alternate approach to carrier tracking applicable to QPSK, consider estimating the carrier frequency offset using the baseband samples after timing recovery (timing recovery with a Gardner TED converges even under relatively large carrier offset frequencies). With that, the correct sample locations for each symbol (just prior to decision) are used with a cross product frequency discriminator to estimate carrier offset as follows:

$$Err = |I_2| |Q_1|-|I_1| |Q_2|$$

Where the absolute values convert all samples to the first quadrant (so this will work for offset frequencies such that the phase doesn't rotate more than $pi/4$ radians between samples, or up to 1/8th the symbol rate.

$Err$ is proportional to frequency offset as the result above is the imaginary portion of the complex conjugate product of the two samples, which is proportional to the phase between the samples. A change in phase versus a change in time (the symbol duration) is proportional to frequency, since $f=d\phi/dt$.

Thus the cross product frequency discriminator is the cross product phase detector for the phase between two samples at two different times.


How is the formula of the cutoff frequency of a 2nd order Butterworth low-pass filter derived? According to articles and notes, the cutoff frequency of a 2nd order Butterworth lowpass filter:

is expressed by the formula:
$\omega _c = {1\over\sqrt{R_1R_2C_1C_2}}$

and since $\omega _c = 2 \pi f_c$, alternatively

$f_c={1\over2\pi\sqrt{R_1R_2C_1C_2}}$

Where does this formula come from? I'm asked to write the equations and formulas which eventually lead to the formula of the cutoff frequency, but I've been very confused and tutorials and notes don't help at all. Thanks in advance.

First you need to find (calculate or look up) the transfer function of that circuit. It has the form

$$H(s)=\frac{a}{s^2+bs+c}\tag{1}$$

where the constants $a$, $b$, and $c$ depend on the values of the resistors and capacitors.

The $3$ dB cut-off frequency $\omega_c$ can be found by solving

$$\big|H(j\omega_c)\big|^2=\frac12\big|H(0)\big|^2=\frac12\left(\frac{a}{c}\right)^2\tag{2}$$

(because $H(0)=\frac{a}{c}$).

Solving $(2)$ involves solving a quadratic equation in $\omega_c^2$.


Power or Energy Signal of x(t) = u(t) Is $x(t)=t \cdot u(t)$ a power signal or an energy signal?
Please "show" me why - preferably using equations.

To calculate the energy of a continuous signal you use the equation:  $$E_\infty = \int\limits_{-\infty}^{\infty} |x(t)|^2 dt$$   and for the power:   $$P_\infty = \lim\limits_{T\to\infty} \frac{1}{2T}\int\limits_{-T}^{T} |x(t)|^2 dt$$  If your signal has finite power ($0<P<\infty$) then it's a power signal. If the signal has finite energy then it's an energy signal (Notice that a signal can be neither type but can't be an energy and a power signal at the same time).


Power of summation of two signals Let $|A|$ denote the power of the analog signal $A$. Is that true to assume $|A+B|=|A|+|B|$, for two arbitrary signals $|A|$ and $|B|$? That is the power of the summation of two signals is -at least roughly- equal to the summation of the powers of the two signals. What if $A$ and $B$ are random signals?

For the amplitude, the inequation $|a+b|≤|a|+|b|$ is satisfied.

However, since the power includes a square calculation on the amplitude, you can find all the three ($=$,$\leq$,and $\geq$) relations.

The difference is that amplitude is a vector, while power is a scalar.


Does this Signal Smoothing algorithm have a name? We are reverse engineering some 20 year old software.  Original developers were laid off years ago, and cannot be found.

In this code, there is a signal that is getting smoothed as follows:

$P_{new} = P_{old} + P_{\Delta} \frac{H + (I-H) P_{\Delta}^2}{J^2 +
    P_{\Delta}^2}$
$P_\Delta = P_{newest Observation} - P_{old}$
$H$ is commented in the code as, Quiescent Smoothing.  It is a hardcoded constant, $H > 0$.
$I$ is commented in the code as Step Smoothing.  It is a hardcoded constant.  $I > 0$.
$J$ is commented as Noise Threshold.  It is a hardcoded constant.  $J > 0$.

I am interested in the name of this algorithm.

I was asked, as I've done a lot of work in motion tracking, but this is clearly not a Kalman filter, or any of it's variants.

Not sure if this has a name, but it is a nonlinear low pass filter that uses different smoothing constants depending on the input signal deviation from the filtered output. Small deviations are typically assumed to indicate consistency with the smooth estimation and result in little adaption to the input, while large deviations indicate a relevant state change and allow the system to adapt quickly.

The system difference equation is
$$y[n] = y[n-1] + \left(x[n]-y[n-1]\right) * F\left((x[n]-y[n-1])^2\right)$$

If you assume that the feedback term is constant and $F=1$, the system passes the input to the output unmodified. Smaller $F>0$ will create a first order lowpass filter and the output signal will be smoothed correspondingly.

If the nonlinear term is well behaved and can be analysed like a time-dependent cutoff parameter. The interesting behaviour of your $$F(d)=\frac{H+(I-H)d}{J^2+d}$$
is that near $d=0$ and for $d\to\infty$. For the first case you get $F(0)=H/J^2$ and for the second case $F(\infty)=I-H$. In between $F$ is monotonous and approximately constant for $d<<J^2$.

That means for small $d$ compared to $J^2$ the filter rejects assumed additive noise with approximately constant feedback and therefore behaves like a traditional low pass filter with feedback $H/J^2$. As $d$ increases, the filter behaviour switches to pass-through for sudden changes and reaches its asymptotic behaviour as traditional low pass filter at a feedback of $I-H$.

Arguably a different set of parameters would make this filter structure clearer. Substituting $H/J^2$ and $I-H$ with frequency related feedback constants would probably help.

Edit: I've used similar filters for smoothing high resolution potentiometer inputs in MIDI controllers.


Pitch detection, YIN, pYIN So here is the "seminal" YIN paper:

  Cheveigne A, Kawahara H. - YIN, a fundamental frequency estimator
  for speech and music

and the new, improved probabilistic YIN:

  Mauch M, Dixon S. - PYIN: A fundamental frequency estimator
  using probabilistic threshold distributions

Now, I am already at a point where I am unimpressed with the only new idea in YIN which is this "cumulative mean normalized difference function" (CMNDF) which seems to have, for its only purpose, to avoid picking the minimum of the "average squared difference function" (ASDF) that always exists at zero lag.  Big deeel.  There are other, better ways of avoiding that obvious error.

What I am trying to figure out with this pYIN is the pitch tracking and pitch candidate evaluation model.  I understand Markov processes and I think I understand what the Hidden Markov Model (HMM) is.  And I believe that they are trying to connect pitch candidates of the previous frame to pitch candidates of the current frame.  I have done this with my pitch detectors (but not so much using an HMM).  We do this to avoid glitches and bleeps arising from an octave error when either the presence of a subharmonic causes the period to appear to be twice as long as the desired period, or weak odd-harmonics causes the period to appear to be half as long as the desired period.  I am quite familiar with this problem and have attacked it with my own bonehead technique which seems to work pretty well, is simple, and makes very few assumptions about the model of the quasi-periodic signal of which we wanna know the period.

I haven't groked from the pYIN paper how a finite (hopefully small) number of pitch candidates are chosen in the current frame and how they are related to pitch candidates of the previous frame.

And then I wanna grok how they eventually score each pitch candidate for the final decision.

In my pitch detectors, I use a variant of autocorrelation computed directly from the ASDF function, and maxima of the autocorrection always correspond to minima of the ASDF.  Avoiding the peak at zero lag is trivial if you guarantee that the input to the pitch detector is DC free, but then of each peak that exceeds a very gross threshold, I handicap peaks at larger lags only slightly so that I do not spuriously pick the lag at two periods simply because it looks slightly better than the lag at one period.

Not shown in that SE answer of mine is how each peak that I consider a candidate to be a pitch candidate (these are peaks that exceed a gross threshold), determine the exact location (to fractional sample precision) and height of the peak using parabolic interpolation, and if that peak is one of the, say, five highest peaks, I call it a "pitch candidate" and associate that candidate with a candidate in the previous frame closest to the current candidate in log frequency (which is pitch).  Then a score is made to each current candidate that is a heuristic function of the current peak height and the previous frame assuming this current peak is tied one of the previous frame.  And for "stickiness", whatever candidate wins in the previous frame, if there is a candidate in the current frame that is tied to the previous frame winner, there is an extra boost to that candidate score.  A slight preference is give to the candidate that I had previously chosen.

But, what is Mauch and Dixon doing?  If there was some code outlining this, I would be interested.  I want to ditch the YIN thing but consider their candidate scoring method.  But I cannot decode it from the paper.

If there was some code outlining this, I would be interested

I don't have the expertise to describe pYin precisely, but I can point you to the original code for that paper https://code.soundsoftware.ac.uk/projects/pyin/repository
There's also a more recent implementation which made it into the librosa library since this question was asked. The code and pull request discussing the implementation are here: https://github.com/librosa/librosa/pull/1063/files#diff-f4245b657ccbb33a7bf56a5f7754585dfbbb80d2e545e840f678aceb548aa206R735


How exactly does sound get encoded by a cochlear implant processor? I've been wanting to know more about how sound is processed by a cochlear implant processor, how they are programmed, and how they stimulate the magnets and electrodes inside.

How do you take pure tones of varying frequencies, noise signals, and complex speech sounds and turn them into a set of organised electrical signals?

Also, what is meant when people say that cochlear implants can have one channel, two channels, etc? I'm only familiar with mono and stereo channel, but recently I heard that channels refer to the splitting of a single signal into two separate paths. If the signals sound exactly the same in both channels, the result is a mono perception. If there is a difference in signal, the result would be a stereo image.

With cochlear implants, they can have up to eight channels or more. How do these work, and just out of curiosity, what sampling rate and bit depth, as well as bitrate is the sound presented? Is it in Pulse code modulation?

I guess the thing to think about is whether the processor, which is digital converts the analogue sound signals of the microphone into a set of signals that somehow need to be reconverted to another form of Analogue that will stimulate an electromagnet that will stimulate a certain number of electrodes for the brain to perceive sound.

There's an excellent non-specialist (science for the public) book on the history of cochlear implant technology:

"I Can Hear You Whisper: An Intimate Journey through the Science of Sound and Language" by Lydia Denworth

which provides a non-technical description of the use of channels in the encoding.


Relation between Normal Phase Shift of a wave and Phase Modulation I am little confused with the Phase Modulation and the phase of a sine wave. I get the phase modulated wave from the google images as below:

Phase modulation has the form:
$$PM(t)=A_c \cos\big(2\pi f_c t+(xA_m)\big)$$
Considering the message signal as
$$ m(t)=A_m\cos\big(2\pi f t\big)\text,$$
and $x$ as the phase modulation index.
Basically, I see $\cos(\theta)$ denotes a wave with some frequency. If we consider in this form - $\cos(2\pi f_1 t)$, it is a cosine wave with frequency $f_1$. And if we add a constant time factor with the unit of seconds, we get something like: $\cos(2\pi f_2 t+t_1)$. I am able to see this as a right phase shifted cosine wave with time $t_1$.
But during phase modulation, we actually change the phase of the wave as per this form:
$$PM(t)=A_c\cos\big(2\pi f_c t+(x A_m)\big)$$

Now what I am confused is that the term $(x A_m)$ has the unit of radians or seconds?
Based on the figure, I am able to assume $x$ has the unit of radians/Volt making the unit of $(x A_m)$ as radians, so that the final phase modulated has an abstract form $\cos(\theta_1+\theta_2)$ which could be some other frequency of the original carrier wave. I am just guessing this but I am not sure. If you could please let me know if this assumption is valid, it would be helpful.

If it is valid, then why do we call it as a phase modulation?
For a phase modulation, I see that the modulated waveform should be something like - $PM(t)=A_c\cos(2\pi f_c t+(x A_m))$ where $x$ should be in the unit of (seconds/volt). Now the phase of the carrier wave changes but I am not able to imagine the waveform.

This does not exist but isn't this the pure phase modulation?

Now, the next hardest part which I am concerned is that, I see that phase modulation is an analog modulation scheme. Let's assume that at a time instant $t_1$, the Phase modulation system calculates the amplitude of the message signal and changes the frequency of the carrier wave proportionally. At time instant $t_2$, this frequency changes to the next value depending up on the message signal.

Now, how fast does the system notice the phase changes of the message signal?
Is that given by the carrier frequency or something else? (Since it is looking the amplitude of the message signal at a specific time instant, can we see it like -the message signal is somehow sampled, its amplitude is quantized to some constant (having the unit of radians) and then added to the angle of the carrier wave?

The reason why I am curious about this is that we can see the phase modulation as the frequency of the carrier wave is changed but that frequency should remain hold at least for one cycle to make it significant.

When adding a constant time shift, the equation is then $\cos(2\pi f (t + t_1))$ (note the parentheses around the time values).
$(xA_m)$ should really be $xm(t)$ in your phase modulation equation - it is a time varying quantity.
$xm(t)$ has units of radians.
It is 'phase modulation' because what is changing directly with the message signal is the phase of the transmitted signal. With amplitude modulation, the amplitude envelope of the transmitted signal carries the message; with frequency modulation, the change in frequency of the transmitted signal (as compared to the carrier) holds the message; and in phase modulation it is the phase of the carrier signal that holds the message. FM and PM tend to look similar because they both affect the argument of the sine wave, but they do so in two different ways.
It is not necessary to have several cycles of the carrier at any given phase shift in order to make that shift significant. Since we know the frequency of the carrier signal, that can be removed from the received signal to recover the message signal. Don't consider it from the perspective of a single time instant (as you point out, this is an analogue scheme, not digital), but rather consider operating on the signal as a signal.
"Now, how fast does the system notice the phase changes of the message signal?"

Continuously. That's how analogue systems work.


Why Does the Median Filter Minimize the Absolute Value Error $L_1$ Cost Function? I can easily prove that the mean filter minimizes the square error $L_2$ cost function using simple calculus.
However, how do you prove that the median filter is optimal with respect the absolute error $L_1$ norm?

Given a set of values $ {\left\{ {s}_{i} \right\}}_{i = 1}^{N} $, we're basically after:
$$ \arg \min_{x} \sum_{i = 1}^{N} \left| {s}_{i} - x \right| $$
One should notice that $ \frac{\mathrm{d} \left | x \right | }{\mathrm{d} x} = \operatorname{sign} \left( x \right) $ (Being more rigorous would say it is a Sub Gradient of the non smooth $ {L}_{1} $ Norm function).
Hence, deriving the sum above yields $ \sum_{i = 1}^{N} \operatorname{sign} \left( {s}_{i} - x \right) $.
This equals to zero only when the number of positive items equals the number of negative which happens when $ x = \operatorname{median} \left\{ {s}_{1}, {s}_{2}, \cdots, {s}_{N} \right\} $.
Remarks

One should notice that the median of a discrete group is not uniquely defined.
The median is not necessarily an item within the group.
Not every set can bring the Sub Gradient to vanish. Yet employing the Sub Gradient Method is guaranteed to converge to median.
It is not the optimal way to calculate the Median. It is given to give intuition about what's the median.


Are 2nd Order Edge Detectors More Susceptible to Noise? I am reading about edge detection and I read that 2nd order detectors are more susceptible to noise. Is there a mathematical proof to this ?

If you assume the Edge Detection is SNR driven operation, one could find a Mathematical justification for this.

First, the variance of Additive White Noise with Variance $ {\sigma}_{n}^{2} $ at the output of a Linear System given by $ g $ is $ {\sigma}_{n}^{2} {\left\| g \right\|}_{2}^{2} $.

Let's look on an Image filtered by a derivative approximating filter $ g = \left[ -1, 1 \right] $.
This increase the Noise Variance by a factor of 2.
For 2nd order the filter is $ {g}_{2} = g \ast g = \left[ 1, -2, 1 \right] $ hence the factor is 6 (It comes from the reason the second time we filter the noise is correlated).

Now, the question does the signal (Edge) is amplified by the same factor?
Moreover, on areas where there is no edge, the probability to detect noise as edge increases.
To see intuitively, just operate 1st order method on the gradient of the image and see how much False edges you have now.


Cropping binary images based on pixel counts in python I have a numpy binary masked image that looks like the following:

I want to automatically crop this image like the following:

Can I achieve this using some sort of image processing rather than building a bounding box regressor?
The method should also work if the image is horizontally flipped.
I am currently using th below crop function to crop the image if any row or column have all 0 pixels:
def crop_image(img,tol=0):
    # img is 2D image data
    # tol  is tolerance
    mask = img>tol
    return img[np.ix_(mask.any(1),mask.any(0))]

This crops out the upper black space from the image but does not work for the lower part for obvious reasons.
My attempt was to try to threshold the image based on pixel counts. Below is the code:
counts = np.sum(image==1,axis=0)
plt.plot(counts)

If we can threshold the image as below we might also consider the task partially successful as we would get something like the following:

Thanks alot. Any help will be highly appreciated.

One approach that is always available for binary images are morphological operations.

The critical piece of code is
# from skimage import morphology
# import numpy as np

# filter elements smaller than 1/6-th of the total width
selem = morphology.rectangle(1, img_binary.shape[1] // 6)
img_opened = morphology.opening(img_binary, selem=selem)

hist = np.sum(img_opened, axis=1)
non_zero = np.argwhere(hist > 0)
top_bound = int(np.round(np.min(non_zero) * (1 - 0.05))) # 5% tolerance
bottom_bound = int(np.round(np.max(non_zero) * (1 + 0.05))) # 5% tolerance

Here is the notebook I used to produce the image above:
https://gist.github.com/FirefoxMetzger/5d0d5e207bf8182f61374ac9ad80f10f


Extraction of non-sinusiodal repetition rates I have an auto-correlation function that was generated from a signal, and I am trying to extract its 'repetition rate' in order to calculate the dominant frequency of the pulse, but I am not exactly clear how to do this.

Here are two cases, labelled 'good' and 'bad' to mark best/worst case scenarios. What methodologies exist that I can use here? Since the repetitions are not sinusoidal, the spectrum does not seem to yield good information.

Here is the original 'good' case:

followed by its auto-correlation.

Similarly, here is the original 'bad' case:

again followed by its auto-correlation:

EDIT: Raw data files as requested:

Raw1

Raw2

EDIT On Feedback:

Spectral analysis was suggested, however I did not have too much luck with it as the spectrum does not look too sharp around the region where you would expect it to. I think this may to due to the fact that the signal is not a repetitive sinusoid and so projecting it on sinusoidal bases does not do too well.

You are right that the repetition is around 650 by how exactly do I compute that automatically? Seems like a peak-picking problem to me? Or is there some other methods that can be used?

Yes, it's just peak-picking.  Your period is the x value of the first strong peak:

Your peaks are all similar in height, probably because you're doing the autocorrelation using the FFT?  So it's a circular autocorrelation.  You can

make it non-circular by zero-padding before doing the FFT, or
could skew the plot you currently have by adding a linear function that emphasizes the peaks close to 0 and reduce the height of the farther ones.

It depends on what you're specifically looking for.  In either case you then have to pick the highest peak that's not at 0 lag.

I adapted my script from here and it (surprisingly) worked without any tweaking:

raw1 x = 709.23

raw2 x = 710.77

from pylab import *

# Import matlab files
import scipy.io
#sig = scipy.io.loadmat('raw1.mat')['raw1'][0]
sig = scipy.io.loadmat('raw2.mat')['raw2'][0]

# Calculate autocorrelation (same thing as convolution, but with
# one input reversed in time)
corr = fftconvolve(sig, sig[::-1], mode='full')

# throw away the negative lags, so x = 0 means 0 lag
corr = corr[len(corr)/2:]

# Find the first low point from derivative
d = diff(corr)
start = find(d > 0)[0]

# Starting at the first low point, find the highest peak
peak = argmax(corr[start:]) + start

# Fit parabola to estimate more precise location of peak
px, py = parabolic(corr, peak)

plot(corr)
axvline(px, color='r')

print px

The variable peak is the integer x value at the peak (709).  Then the parabolic function fits a parabola to the peak and tries to get a more precise estimate (709.232...), which may or may not be relevant to what you're doing.

The python fftconvolve() function does zero-padding internally to produce a linear convolution, like convolve() does


How can understand periodicity of a Signal from frequency domain representation? Is it possible to say a signal is periodic from its frequency domain representation? A periodic signal is sum of its sinus and cosinus. Frequency translation of sinus and cosinus functions are symmetrical impulses. Is it correct to say "a signal is periodic if its frequency domain representation includes symmetrical inpulses"?

A signal is periodic if and only if its frequency spectrum is

discrete (contains only pulses covering a single locus with a finite
area).
all of those discrete pulses occur at an integral multiple of one
frequency (the periodicity).

If you have, for example, $\cos(\pi t)+ \cos(t)$, the resulting spectrum is discrete but the function is not periodic.


Indetermination in a complex fourier series I determined the complex Fourier series of a sinusoidal signal and arrived at the following expression:
$$\sum_{n=\infty}^{\infty} \left[\frac{4e^{-j \frac{\pi}{2}n}}{\pi(1-n^2)}(e^{-j\pi n}+1)\right]e^{-j \frac{\pi}{2}tn} $$
However, we have a division by zero at $n = 1$ and $n = -1 $. Is there any way to determine $D_{-1}$ and $D_1$?

HINT:
You need to go back to the original definition of the Fourier coefficients via the integral, and treat $|n|=1$ as a special case. Your mistake happened there already (you divided by zero).


Need samples of big endian encoded tifs Folks I've written a small utility that extracts some (and only some) header information from Tifs.

To test it fully I need some Big Endian Tifs, but I cant find any, and I can't find anything to generate them.

Does anyone know of a public repository for tifs with big endian encoding? Or of something that can generate some?

Thanks, apologies if this is off topic.

The current ImageMagick can make them.  Given a little-endian TIFF,

convert le.tiff -endian msb -define tiff:endian=msb be.tiff

and vice versa:

convert be.tiff -endian lsb -define tiff:endian=lsb le.tiff

Run

identify -verbose be.tiff | grep ndian

to see what happened.


Controlling one signal with another one (Rocket engine effect - Subtractive synthesis of sound) I am programming some sound effects in Java and exporting them into .wav files. Currently, I am trying to program a rocket engine sound effect. I want to do it in the following way:

The sound of a rocket engine may be synthesized with a red noise generator controlled by a second red noise generator. The parameter of the first generator modified by the second one is the number of interpolated samples, influencing the spectral content of the generated noise. In order to simulate changes of sound intensity (e.g. during launch) the envelope generator should be used.

I am wondering how can it be done, e.g. what does it mean that one signal controls another one. Probably this part is explaining it, but I am not sure what to do now:

  The parameter of the first generator modified by the second one is the
  number of interpolated samples, influencing the spectral content of
  the generated noise.

Is it about this parameter describing how many values from white noise are taken and linearly interpolated while creating the red noise? (see my simple drawing explaining this process below)

I have a red noise generator, which returns an array of doubles with values between -1 and 1 (it is generated from the white noise as described). What am I supposed to do now? How can I control the second red noise? I guess that it does not mean that I should control the amplitude of the second signal. Does it?  Schema of steps required to obtain the rocket engine sound effect is attached below.

I have already solved this problem:
I have a white noise generated as an array of doubles:

public static double[] whiteNoise(int length) {
    double[] out = new double[length];
    for (int i = 0; i < length; i++) {
        out[i] = (Math.random() * 2 - 1.0)/2.0;
    }
    return out;
}

I obtain the red noise from the white noise in the following way:

public static double[] redNoise(int length, int m) {
    double out[] = whiteNoise(length);

    for (int i=0; i<length/m-1; i++) {
        int k = i*m;
        int l = ((i+1)*m <= length-1) ? (i+1)*m : length-1;
        double d = (out[l] - out[k])/((double)(m));
        for (int j=0; j<m; j++) {
            out[k+j] = out[k] + d*j;
        }
    }
    return out;
}

Where the m parameter is the number of interpolated samples (we take the first sample, and the m+1 sample, and interpolate values of all samples between these two).

And now the SHOW BEGINS

I generate the red noise with the interpolation number equals to let's say 10000. According to the value of consecutive samples, I alter this red noise with different values obtained from these samples taken as the interpolation number parameter. In this way a rocket sound effect is obtained.

int interpolatedframes = 10000;
    double[] rednoise1 = SigGen.redNoise(duration, interpolatedframes);
    while(i<rednoise1.length-1)
    {
        selectedNumber = (int)Math.floor(Math.abs((100*(rednoise1[i]+0.5)+100)));
        if(i+selectedNumber < rednoise1.length-1) {
            rednoise1 = SigGen.alterRedNoise(rednoise1, i, selectedNumber);
        }
        else
        {
            break;
        }
        i+=selectedNumber;
    }

With the following code the array with the noise is altered in a way, that it interpolates "interpolationNumber" of samples beginning from the startingIndex.

public static double[] alterRedNoise(double[] input, int startIndex, int interpolationNumber)
{
        int k = startIndex;
        int l = ((startIndex + interpolationNumber)<input.length-1?(startIndex+interpolationNumber):input.length-1);
        double d = (input[l] - input[k])/((double)(interpolationNumber));
        for (int j=0; j<interpolationNumber; j++) {
            input[k+j] = input[k] + d*j;
        }
    return input;
}

Sorry for a bit messy code, but I just wanted to publish it ASAP.

#EDIT 1

The sound which is supposed to be generated is the sound of a rocket being launched. Take a look at this movie (maybe from 1:45): https://www.youtube.com/watch?v=OnoNITE-CLc

To have the effect of the rocket flying over, an envelope generator might be used. It can be implemented as a simple multiplication of the samples obtained from the described method in such way, that it will be gradually suppressed.


How to extract white objects in an image(attached)? The program should then display each N as a separate figure. The program should take the given image as input and the end result should be each white object in a separate figure window.

It's all about detecting intensities of pixels. It'll be better if you use MATLAB. Such binary images will have just 2 intensity levels 0(black) and 1(white).

No one is going to do it for you unless you are ready to pay. Buy and read a copy of Digital Image Processing by Gonzalez


Why does signaling overhead for time synchronization scale up with the number of transmitter nodes in a multiple access system? I was recently reading about random multiple access methods without time synchronisation (of the pure Aloha type). In this paper it is stated that,

  The need for transmitter synchronization is a major drawback for large networks as the signaling overhead scales up with the number of transmitters independently from their traffic activity factor.

I understand that oscillators tend to drift in time and present instabilities depending on the technology, operating conditions, etc., and therefore a reference signal must be given by the central node with a periodicity that depends on the desired accuracy and the transmitter's clock drift.

I would like to know why there is also a dependency on the network size.

EDIT: The only scenario I can think of where signalling scales up with the network size, is a distributed network (but in that context, we usually talk about logical clocks and not about hardware clocks/oscillators). Note also that the network described in the paper has a star topology, so a periodic broadcast message from the central node would do, independently on the number of users in the network.

Let's look at it a bit more in context:

  Moreover, E-SSA (as SSA) has another major advantage over SA/DSA/CRDSA
  i.e. its truly asynchronous RA nature eliminates the need to maintain
  accurate slot synchronization among all transmitters unlike the case
  in SA/DSA/CRDSA.

  The need for transmitter synchronization is a major
  drawback for large networks as the signaling overhead scales up with
  the number of transmitters independently from their traffic activity
  factor.

  The capacity of asynchronous collision channel without
  feedback was investigated by Massey in his seminal paper together with
  protocols sequences for achieving the capacity boundaries [17].

If I understand well, this is not about frequency locking, but about time slot synchronization: the synchronization occurs on a time scale that is orders of magnitude higher than the oscillator period.

In time division multiplexing, there is a guard period at the end of each time slot where nobody should transmit, both to prevent interference and to allow synchronization.
Because this guard interval is present at each time slot, the sun of guard interval times can be considered to be an overhead that varies with the number of time slots (i.e. for the max. number of users the system is designed for).


Interpolating/Decimating CIC Filter Group Delay I have an Interpolating Cascaded Integrator Comb (CIC) filter designed as so:

Rate Change = 100
Number of Stages = 4
Differential Delay = 1

I'd like to calculate the group delay of this filter but I am not sure how.  This question on DSP Related is helpful, but I am not sure my filter will have 0 delay.

Also, Richard Lyons explanation is helpful, but he doesn't directly deal with group delay.

I also have a Decimating CIC filter designed as so:

Rate Change = 25
Number of Stages = 4
Differential Delay = 1

How would I calculate its group delay?

Consider an $D$-tap FIR filter with liner phase, the group delay (measured in samples) is $$g=\frac{D-1}{2}\tag{1}$$
and therefore, if it is measured in seconds it will be $$g=T_s\frac{D-1}{2}\tag{2}$$
where $T_s=1/F_s$.
The CIC filter which is also denoted as recursive running sum filter is indeed a special implementation of a moving-average filter. The moving average filter is
$$y[n]=\frac{1}{D}[1+z^{-1}+z^{-2}+\cdots z^{-D+1}]=\frac{1}{D}\sum_{n=0}^{D-1}z^{-n}=\frac{1}{D}\frac{1-z^{-D}}{1-z^{-1}}$$
which can be implemented by a $D$-tap FIR filter. If you compare carefully, the above response is identical to the frequency response of each stage of a CIC filter.
Therefore, we can use $(1)$ and $(2)$ for calculation of the group delay in each stage.
Comparing the response given for moving average with that of the CIC filter, assuming $R$ is rate change and $M$ is the number of samples per stages, then $$D=RM$$ So for $N$ stages we have $$g_{\text{cic}}=N\left(\frac{RM-1}{2}\right)$$ expressed in samples and
$$g_{\text{cic}}=NT_s\left(\frac{RM-1}{2}\right)$$ in seconds.


What kind of filter is this? Is this filter a BPF?
$$\dfrac{z}{z-a}$$

where $a$ is some complex number?

If we put a pole somewhere on the unit circle it will emphasise a certain frequency, is that right?

Yes, if you put a pole approaching the unit circle (but not on the unit circle, as this will create an unstable filter) , you will have a bandpass. The center frequency of this bandpass filter depends on how many radians along the unit circle where you have placed the pole (the pole's 'degree').

Moving counter-clockwise from zero will give you an increasing center frequency of the bandpass filter.

Once you start moving towards zero on the Z-plane and away from the unit circle, you will decrease the sharpness Q and consequently increase the bandwidth.


Confusion in Channel Encoding and Convolutional Encoding If I'm making NRZ signal by first encoding 1 input bit into 10 bits multiplied by -V/+V before sending it into a channel, is it called channel encoding?
If I then modulated it using GMSK technique, can it be considered Convolutional Encoding?

Like other comments, I don't really understand your question. What I am trying to do is to list basic things so that other ones can suggest edit because I find it is much easier to write in answer part :). Please let me know if it makes you less confused.

The general communication system is

At transmitter, before modulator is digital, modulator is the block to convert from digital world to analog world because physical signal is transmitted in analog (at least at our normal scale). Modulator converts bit sequences to voltage (NRZ), amplitude and phase of sinusoid (QAM), frequency (FSK), phase of of sinusoid (MSK, GMSK), ...

Signal is transmitted in physical channel always suffers degradation (noise, fading, ...) which is random. To combat this negative effect, channel coder adds redundant information to help us to recover the correct signal. For example adding 9 supplementary bits instead of sending only one, there is the chance that more than 5 bits are correctly decoded and with majority rule, the transmitted data is recovered. The operation is called channel coding. There are a lot of channel codes : turbo code, LDPC, convolutional code, ...

Thus, channel coding and modulation are two different steps that have different roles. For more details, you need to read the recommended resources in the comments.

It seems that you understand "convolutional coding" as modulation. Maybe it comes from the fact that signal transmitted over channel is modeled as convolution operation ? People usually use "convolutional code" to talk about a specific family of code in channel coding step.

About the source coder block in the figure, on the contrary of channel coder that adds redundant information to combat random degradation, source coder is here to remove redundant information from source. Huffman code, Shannon code, ... are famous. Applications such as winrar, winzip, ... are well-known. Source coder and channel coder have different roles, but they are related in the sense of "optimality" because they both deals with "information". Thus it is not always correct to say that source coder and channel coder are independent. You can search for "joint source channel coding" if you are interested. The reason why the myth exists and an counter-proof can find in this article Source Channel seperation Theorem. Note that it could be an advanced topic, depending on your background.


What are thoughts on best modulation to use for underwater acoustic communications at low frequencies? I wanted to ping the DSP hivemind for general thoughts on what would be the best modulation type to use for low frequency underwater communications. I have chosen this project since I can learn a lot from it.

Some context:

Low frequency as in < 500 Hz (for both carrier and modulated data obviously)
BPS of say, 200 Hz would be good.
Will certainly have multipath.
Frequencies can be smeared due to doppler by a maximum factor of about 0.3% of the original frequency.

What I have found thus far:

I was thinking about OFDM, but I have learnt that while channel estimation is much easier, it is much more sensitive to doppler effects.
I was also thinking about chirp modulation, has anyone ever done something like that?

What are your thoughts?

Edit: I have attached some of what I believe are 'worst case' scenarios of the (multipath channel, for a bps = 200 Hz). The channel is in terms of number of bits in the time domain so you can more easily see how many bits pass before the next reflection arrives.

Case 1:
Case 2:
Case 3:
Case 4:

Notes:

As we can see, I almost always have a second path of nearly equal magnitude but opposite phase ready to stick it to my main path.
For a packet of 1000 bits @ 200 bps, (5 seconds), I suppose the channel could change significantly... but at the same time, we have total control over the packet length and content.
We may assume that frequency offsets due to doppler is relatively 'well behaved', ie, no sudden 'jerks'. Frequency offsets due to carrier mismatches can also be thought of as well behaved.

That is a nasty signal environment.  I would do a little DSP judo and make the multi-path work for you by using a rake receiver, which should make the multi-path signals increase your SNR rather than diminish it.

Rake receivers have only been used (as far as I know, anyway) in CDMA systems, but that is not because they can only be used in CDMA systems.  The issue is the symbol period.  The following quote is from a paper on rake receivers:

  RAKE receiver attempts to
  collect the time-shifted versions of the original signal by
  providing a separate correlation receiver for each of the
  multipath signals. This can be done due to multipath
  components are practically uncorrelated from another when
  their relative propagation delay exceeds a chip period.

Thus, the reason that they were used in CDMA systems and not, say, GSM, was because the chip speed was so much faster, and thus the multipath signals would not arrive in the same chip/symbol.  That was not true for GSM and other "narrowband" signals.  Even though your bit rate is very low, the multipath signals still seem to have delays greater than one symbol, so this constraint should not be a problem.

The other problem is detecting the multipath.  This could, I would think, be accomplished with known data sequences, such as a preamble, but I confess that I am not an expert on rake receivers.

If the rake receiver does not work out, you could still use a long equalizer to handle the multipath.

Regarding the other elements of the system, if there is no feedback I would do a QPSK signal with FEC (probably turbo codes).  If you do have feedback I would do the same but also dynamically change the modulation type from QPSK up to 16-QAM or something, depending on how well the signal is coming through.

EDIT: After thinking about it a bit I realized why rake receivers and CDMA systems go hand-in-hand.  The problem is that even if you detect a multi-path signal it doesn't do you a lot of good unless it has a positive SINR.  By definition at most one multi-path signal can have a positive SINR, because for all the other multi-path signals the strongest overwhelms them.

That is where despreading comes in.  Once the multi-path signal is despread it will have a positive SINR, assuming the spreading factor is large enough to overcome the initial negative SINR.  Given this I think the right solution is to relax the 500 Hz limit, use a fairly large spreading factor and a rake receiver to combine the various multi-path signals.


Spectrum analysis of a function I asked this on Math SE but received no replies. I hope this is a relevant forum to ask in.

I would like to analyze the spectrum of the following function:
$$f(t)=\cos(t\cdot a(1+b\cos(ct)))$$
with $a,b,c$ some real parameters. As the function is non periodic (perhaps unless $a,b,c$ are fine-tuned), the most natural approach is to Fourier transform it. But doing it by hand is quite horrible, Wolfram Alpha aso struggles, and using Mathematica I get
$$-ia\cdot\cos(\sqrt{2\pi})\cdot \delta'[\omega] - abc\cdot\cos^2(\sqrt{2\pi})\cdot \delta''[\omega]$$
I'm only interested in the real part of the result. I'm now trying to make sense of what I got.

1) Does the Mathematica result seem reasonable to you?

2) Looking at it from a generator function perspective (as in here), it seems there is a sharp peak centered at zero, accompanied by more minor peaks at $\pm\Omega$ for some $\Omega>0$. I find this result surprising, as I expected that the spectrum would be centered at $a$ as the "main" frequency followed by some smaller harmonics or some sort of $sinc$-like function.

3) Following a much simpler argument than the full Fourier transform, I try to think of it as a function that has a frequency that oscillates in "time" ($t$). Then maybe as an analogy I can ask myself which frequencies are most probable. In that case, I can think of the $\cos(ct)$ term as the cosine of a uniform random variable $ct:= Y\sim Uni(-\pi,\pi)$, and so the probability density function is actually U-shaped centered at $a$ and given by a function in the general form of $$f_Y(x)=\frac{1}{\pi\sqrt{1-x^2}}$$(normalized etc.)

But these two analyses are far off from one another, and I would like to understand what is going on.

A counter-argument to my probabilistic interpretation is that the term that determines the intensity of a frequency component is the coefficients of the Fourier decomposition (if it were periodic), and here this is not a linear combination of frequencies with different probabilities, but something else. Still, I need some help explaining the results. Thanks!

The phase term of your outermost cosine function is time-varying: $\phi(t):=t a (1+b \cos(ct))$. The instantaneous frequency is given by $\frac{\partial \phi}{\partial t} = a+ab \cos(ct) - abct \sin(ct).$ In particular, in a tiny neighborhood around $t=0$ we can argue that $f(t)$ has only one frequency component $a+ab$. At large values of $t$, the $abct\sin(ct)$ term takes control.

Notice that the frequency varies with time. This makes it difficult to intuitively interpret the Fourier transform which assumes that all frequencies exist for all times. For such signals, it may be better to use time-frequency analysis by looking at the spectra of windowed segments of $f(t)$ around different values of $t$.


What are Chunks when recording a voice signal I am having trouble understanding the example program for audio recording ,
it comes with PyAudio module.

I am using an external Audio Card to record the sound , obviously, the voice should be converted into the digital form to be presented to the computer. But i do not understand how this code converts the input Analog Voice into Digital Form and then saving it.

Here is the code,

import pyaudio
import wave
import sys

CHUNK = 1024
# What is CHUNKS here ?
FORMAT = pyaudio.paInt16
CHANNELS = 2
RATE = 44100
RECORD_SECONDS = 5
WAVE_OUTPUT_FILENAME = "test.wav"

p = pyaudio.PyAudio()

stream = p.open(format=FORMAT,
                channels=CHANNELS,
                rate=RATE,
                input=True,
                frames_per_buffer=CHUNK)
print("start....")

frames = []

for i in range(0, int(RATE / CHUNK * RECORD_SECONDS)):
    data = stream.read(CHUNK)
    frames.append(data)

print("done...")

stream.stop_stream()
stream.close()
p.terminate()

wf = wave.open(WAVE_OUTPUT_FILENAME, 'wb')
wf.setnchannels(CHANNELS)
wf.setsampwidth(p.get_sample_size(FORMAT))
wf.setframerate(RATE)
wf.writeframes(b''.join(frames))
wf.close()

When the code is executed it record the sounds for 5seconds using an external audio card with a mic.

What I understood from the code is,

The Rate is 44Kbps, which means that 44 kilo bytes are generated from the Analog Audio signals per second.

Then an empty frame is generated frame[] where all the recorded data will be stored.

The next thing is for loop,

for i in range(0, int(RATE / CHUNK * RECORD_SECONDS)):
    data = stream.read(CHUNK)

the answer to RATE/CHUNK*RECORD_SECONDS is 215.33 which means there will be 215 iterations.

I don't understand why the RATE is divided by CHUNKS and then multiplied by RECORD_SECONDS for the recording of 5 seconds. (RECORD_SECONDS = 5)

Also, What exactly is CHUNK ?

So essentially, @porten pointed out what chunks of sound are.

Let's look at your code:

CHUNK = 1024
# What is CHUNKS here ?

The chunk is like a buffer, so therefore each buffer will contain 1024 samples, which you can then either keep or throw away. We use CHUNKS of data, instead of a continuous amount of audio because of processing power. (Let's assume this was a continuous flow of data being received from a Microphone and is being recorded and saved) then it would just eat up the processor, thus causing potential crashes. In terms of Raspberry Pi / Arduino development (Where the RAM is very small) CHUNKING the data like this makes the stream flow easier and thus prevents memory leaks.

Another reason:

Let's assume that you wanted to implement an algorithm for determining whether something is speech, or, is just noise. How would/could you do this using a constant stream flow of sound data? It would be very difficult. Therefore, by storing this into an array (or list in your case) you can perform analysis on this data, for example RMS. Then you could have some threshold to determine if you want to keep the data, or the data is no good. E.g.

CHUNK = 1024

var input[CHUNK] = readAudio(CHUNK); // reads in 1024 of data

// compute RMS
a = RMS(input);
if(a >= THRES_VALUE)
{
    actual += input;
}

So in this example above, we record 1024 bytes of data, we then compute RMS and check to see whether or not the sound is speech or just white-noise.

To confirm as well, the chunks contain the data.


Quadrature subsampling result not as expected Let's suppose I have an ADC running at 81.824MHz. I am quadrature subsampling a signal that is at 143.5MHz using a sampling rate of 81.824MHz/4 = 20.456MHz = Fs. I expect the signal to appear at baseband with a small offset of 143.5MHz-7*Fs = 308KHz. But it is showing up at 320KHz when I run the code below. Why is that?

% Sampling frequency
fs = 20.456;

% ADC frequency
fadc = 4*fs;

% ADC period
T = 1/fadc;

% Some ADC samples
n = 0:T:100;

% Received frequency
f_l1 = 143.5;
s_l1 = cos(2*pi*f_l1*n);

% Replica frequency for mixing.
f_local = 7*fs;
s_local = cos(2*pi*f_local*n);

% Mixing
sif = s_l1 .* s_local;

% Find frequency content
siffft = (fft(sif));

% Plot results
nyquist = 1;
freq_steps = (1:length(siffft))/(length(siffft))*nyquist*fadc;
figure(3);
plot(freq_steps, abs(real(siffft))),
grid;

Your calculation is correct; sif has a harmonic at 308 kHz. Given the way you have set up the problem, you can't get that exact value using the FFT because you don't have a frequency bin at that frequency. According to my calculations, the nearest bin is at 305 kHz, and I do see a peak at that frequency. Since there is no bin at the signal's frequency, the energy will leak: bins in the neighborhood of the peak will also get some of the energy of the 308 kHz signal.

Leakage is very hard to avoid. Assuming a fixed sampling frequency, you can calculate the number of samples you need to have a bin at exactly 308 kHz. However, it's very likely that you will then have a non-integer number of cycles of sif, so you will have leakage again. This effect can be reduced by using windowing.

There are many related questions in this website (for example, this one). Search for 'fft bin leakage' to find many more.


Linear combination of DT unit impulse Came across this example in class but I'm not sure how the expression $x(n)$ was derived.

The plot that shows the sequence $x[n]$ is also called as the lollipop plot and it shows the individual samples of $x[n]$ with their weights and locations.

Each lollipop represents a discrete-time unit-impulse signal, $\delta[n]$, which is located at the position it's shown and weighted by the value of the signal there. The superposition sum of all those impulses compose the signal $x[n]$.

Now, your signal has 3 lollipops located at $n = \{-1,0,1\}$ and weighted by $\{1,2,3\}$ respectively.

Since an impulse located at the position $n=d$, and weighted by $a$, is given by
$$ a \cdot \delta[n-d]$$

by using this, you can compose the signal $x[n]$ from its lollipops as:

$$ x[n] = 1 \cdot \delta[n-(-1)] + 2 \cdot \delta[n-(0)] + 3 \cdot \delta[n-(1)] $$

$$ x[n] = 1 \cdot \delta[n+1] + 2 \cdot \delta[n] + 3 \cdot \delta[n-1] $$

Note that a generalization of this also yields the sifting-sum decomposition of any signal $x[n]$ as

$$x[n] = \sum_{k=-\infty}^{\infty} x[k] \delta[n-k] $$


Is there a $z$-transform like for variable sampling rate signals? I'm working with signals with variable sampling rate (the time space between samples is not constant). I know the delay between samples but I don't wont to interpolate the signal.

Is there a way to design a filter in $z$ domain and filter a signal with variable sampling rate without interpolation?

Since the powers in the $z$-transform relate to delays or lags, in the philosophy of Laurent series, I would like to mention Puiseux series, a generalization of power series with negative and fractional exponents, or the Hahn (or Hahn–Mal'cev–Neumann) series, with even more flexibility in the  exponents. However, their use in signal processing remains quite shallow to me yet.

I suspect than rational exponents could be still tractable. Indeed, they are used on multiband ($M$-band) filterbanks, where you can spot weirdities like $X(z^{1/2})$ or $X(z^{1/M})$, but most of them are based only on rationals with $M^p$ denominators.
More generic powers (irrational) could be way more involved. The only mentions of Puiseux series I have found on IEEE Xlore are in the domain of control and automation (definitely not my cup of tea).


Sinusoid with increasing frequency How would you describe this signal? It's like a sinusoid but as if its frequency was constantly increasing: could you write down a mathematical description? Thanks. And no, this is not homework...

Frequency is more likely decreasing (period is getting longer). This is a sweep sine signal or as some people used to call it - chirp tone. No point in rewriting the literature, so here are some links with very good explanation for both linear and exponential frequency change.

For linear sweep you get:

$$s(t)=\sin\left(\omega_1t+\frac{\omega_2-\omega_1}{T}\cdot \frac{t^2}{2}\right) $$

where:

$\omega_1 $, $\omega_2$ - initial and final angular frequency respectively

$T$  - sweep duration.

  Meng Q., et al. - Impulse Response Measurement With Sine Sweeps and
  Amplitude Modulation
  Schemes

  Farina A. - Advancements in impulse response measurements by sine
  sweeps


3D world homography matrix calculation I have a matrix in (x,y,z) coordinates and also i have a transformed version like (x_new,y_new,z_new).

I need to calculate a transformation matrix that can transform (x,y,z) to (x_new,y_new,z_new). But i have really no idea is there any kind of technique?

Check out my answer to a somewhat similar problem here.  You will need to generalize

$$ x = A_x m^2 + B_x mn + C_x n^2 + D_x m + E_x n + F_x $$

to

$$ x_{new} = A_x x^2 + B_x y^2 + C_x z^2 + D_x x y + E_x x z + F_x y z + G_x x + H_x y + I_x z + J_x $$

Also the same for $y$ and $z$.  This will lead to needing to solve a 10x10 inverse.  If it is linear, which will be confirmed if the second order coefficients are all zero or near zero, then you can put your results in matrix form like this:

$$
\left[
\begin{array}{c}
x_{new} \\
y_{new} \\
z_{new} \\
\end{array}
\right]
=
\left[
\begin{array}{cccc}
G_x  &  H_x  & I_x & J_x \\
G_y  &  H_y  & I_y & J_y \\
G_z  &  H_z  & I_z & J_z \\
\end{array}
\right]
\cdot
\left[
\begin{array}{c}
x \\
y \\
z \\
1 \\
\end{array}
\right]
$$

If you prefer, this will also work:

$$
\left[
\begin{array}{c}
x_{new} \\
y_{new} \\
z_{new} \\
\end{array}
\right]
=
\left[
\begin{array}{ccc}
G_x  &  H_x  & I_x \\
G_y  &  H_y  & I_y \\
G_z  &  H_z  & I_z \\
\end{array}
\right]
\cdot
\left[
\begin{array}{c}
x \\
y \\
z \\
\end{array}
\right]
+
\left[
\begin{array}{c}
J_x \\
J_y \\
J_z \\
\end{array}
\right]
$$

If it is linear, recalculate with this equation and its $y$ and $z$ counterparts:

$$ x_{new} = G_x x + H_x y + I_x z + J_x $$

Then you will only have to solve a 4x4 inverse and will get a better fit.

In the 10x10 case, you will need at least 10 samples.  In the 4x4 case, at least 4.  The more you have, the better you can tell if your results are good.  This technique produces a best fit solution.

It is also possible to generalize to higher powers, but of course this will lead to much larger inverses needing to be solved and more samples.  You can think of this approach as a generalized multidimensional Taylor series best fit.

Hope this helps,

Ced


Sensor fusion with accelerometer and gyroscope only I need to track an object in 3D-space and I can use only accelerometer and gyroscope. I can't use a magnetometer because the object itself has a strong magnet.

I need to know the orientation of the object (how much gravity there is in $x, y,$ and $z$ components), and how much acceleration there is in $x, y,$ and $z$. I'm not trying to calculate the object's absolute position.
Is it possible to create sensor fusion that has less than $10\%$ error in gravity and acceleration? At this point I'm just trying to understand if this goal is realistic.

I did a quick experiment with my smartphone and found out that the way it calculates gravity (LPF?) results $\approx 40\%$ error when phone is moved.

This was just a test: I'm not making phone app. I'm designing dedicated HW.

I'm just posting here a pic to show how the "gravity sum" looks like after fixing the math as suggested by Maximilian.

Actually, the plots for the components of gravity look good. What you miss is the fact, that your overall gravity (i.e. gravity sum) is the length of the vector, that is spanned by the three components. Hence, you should calculate

$$G=\sqrt{g_x^2+g_y^2+g_z^2},$$

where $g_x,g_y,g_z$ are the three single components measured by the acelerometer. Looking at your curves, I predict the sum will be quite flat.


This shows the subset of the data you included and the absolute value of the DFT (discrete Fourier transform) of a zero-padded version of it.

R Code Below
    #Q31040

data <- read.csv('Q31040/data.csv')

par(mfrow=c(2,1))
plot(data, type='l', col="blue")

zero_padded_data <- c(data$Volts, rep(0,1000))

omega <- seq(0,2*pi - 2*pi/length(zero_padded_data),2*pi/length(zero_padded_data))
fft_data <- abs(fft(zero_padded_data))
plot(omega, fft_data, type='l', col="blue")

peak <- which.max(fft_data)
points(omega[peak], fft_data[peak], col="red", lwd = 5)
title(paste("Peak location is", omega[peak]))


How to represent ratio between harmonics in audio spectrum I was provided a report where the amplitudes of spectral harmonics have been extracted from several audio files. The amplitudes have been extracted using Audacity's spectrum visualizer, where the spectrum is provided in dB. I want to compute the ratio between the amplitude of the fundamental and the amplitude of each harmonic. Does it make sense to perform a ratio in dB or should I use the difference between amplitudes? Thank you.

I think what you are referring to is called the Total Harmonic Distortion (THD). THD is the level of the fundamental divided by the level of the harmonics.  You can express this as either as a ratio (X/Y), a percentage (X/Y * 100) or in dB (20*log10(X/Y)). A percentage is probably what most people are familiar with.


low pass filter with minimal ringing I am trying to filter out some square wave signal to within a limited band (1/4 or 1/8 of the original), I realized that there's a lot of ringing in the wave when I use my filter (elliptical), I also tried Butterworth, and others (given in Matlab fir1, and classic iir filters) but the only filter that seems to give no ringing is Gaussian. So my question is, how should I go about designing a LPF with minimal ringing? (preferred characteristics: minimal pass band ripple, stop band of more than -50dB, relatively fast roll off). Also as I am trying to implement this in a DSP, low order filters such as IIR types are preferred.

Thank you for any help.

Generally the amount of ringing that you get is a function of the steepness of the filter in the frequency domain, regardless of filter type.
At the same steepness an Elliptic will require a lower order but will have pretty much the same ringing as a Butterworth.

Choosing between a linear phase or minimum phase will change the character but not the extent of the ringing.

No matter how you slice it, a square wave contains a lot of high frequencies and once you filter those out, it's not a square wave any more.


Filtering on a time signal I am new to DSP, and am mostly self taught. I have a couple of questions, that may be too basic, but I couldn't really find an answer online.

I have a signal obtained by plotting the motion of a human in 2D from a fixed point of reference (origin).

The data is noisy, so I wanted to use a filter to get a smoother signal.
I want to filter out different motions. The intuition is that different kind of motions would occur at different frequencies.

From what I googled, and read online; it seems I need a bandpass filter like butterworth. However, after reading a few tutorials and book chapters, I have two questions:

My signal is a time-signal (displacement vs time), but the filter takes the input in terms of cut-off frequency. Does this mean I have to convert the signal into frequency domain, or does the filter take care of it? I am bit confused on how this actually works.
On a related note, I don't know the frequency of each motion that I am trying to filter out i.e. the intuition is that the frequency for different motions be different, but how do I determine by cut-off values? Is it usually done by trial and error, or is there a more scientific way to accomplish this?

My signal is a time-signal (displacement vs time), but the filter takes the input in terms of cut-off frequency. Does this mean I have to convert the signal into frequency domain, or does the filter take care of it? I am bit confused on how this actually works.

You design or specify the filter in the frequency domain, but the actual filtering can be directly accomplished in time.  There are FFT based filtering algorithms that are typically more efficient, but there is no requirement to transform your signal to the frequency domain.

On a related note, I don't know the frequency of each motion that I am trying to filter out i.e. the intuition is that the frequency for different motions be different, but how do I determine by cut-off values? Is it usually done by trial and error, or is there a more scientific way to accomplish this?

Anything that doesn't require trial and error is suspect, scientific methodology wise. I suggest that you try to be more specific about what filtering out different motions actually means.  Some kinds of phenomena are a superposition of features which are often separable, of which some can be "filtered out" and others don't lend themselves to that.

You might want to use a tool known as the Short Time Fourier Transform (STFT) and see if you can eye ball any patterns associated with desired/undesired motion.

Finally, try Google Scholar.


The Distribution of Filtered Gaussian White Noise Suppose I have some time series $s(t)$ which contains Gaussian white noise generated by a distribution $N(0,\sigma^2)$

Then I apply a filter to s(t) with a frequency response $H(\omega)$, giving me $s_H(t)$

What distribution does $s_H(t)$ have?

I'm particularly interested in cases where $s_H(t)$ is band limited Gaussian white noise.

If you filter a Gaussian random process with an LTI system, the output will also be Gaussian. You can make intuitive sense of this by considering that a linear combination (which is what filtering does) of jointly Gaussian random variables is a Gaussian random variable.

You can find an in-depth treatment of filtering random processes in this MIT OpenCourseWare document (section 7.4).


Nyquist Frequency and Window Length The Nyquist theorem states that the sampling rate must be twice the highest frequency to be observed. How does the length of the interval play into this relationship?
The main resource that I found here states that

For a waveform to be fully characterized, it must complete at least one full cycle within the duration of the time record.

Is this the only restriction, assuming that my sampling rate exceeds the Nyquist rate? For example, assume that the frequency band I wish to measure is from 8Hz to 13Hz. I have a sampling rate of 2000 Hz, which is way beyond the Nyquist frequency. As long as my time window is at least 1000ms / 8 Hz = 125 ms, should I be able to reconstruct the signal? For context, I am planning on running an FIR filter followed by a Hilbert transform on these windows to extract phase information, so I want the analytic signal to be mathematically sound.

(this should possibly be a comment but I don't have the rep yet, so I'm submitting it as an answer) Restricting your view of the signal to a finite time window will always introduce artifacts, unless the signal itself is zero everywhere outside of that window. Using a "square" window (simple sampling the signal for a finite duration) will have the effect of convolving the signal's true spectrum with a Dirichlet Sinc function (Dirichlet Kernel). See https://en.wikipedia.org/wiki/Window_function#Spectral_analysis


Identify the Type of Image Distortion (On Lena Image) I want to understand the distortions in each of the lena images shown in the image below.

I see that the first one is a blur. How can undo the blur. The second one is colored salt and pepper noise. This is different from the gray salt and pepper. How can remove that from the image. And for the last one, I don't know that kind of distortion is in that image.

Blur
If you want to reverse Blur applied on an image (Using Convolution, namely Linear Spatial Invariant Blur) you should use Deconvolution which, as name suggests, the inverse operation of Convolution.
Have a look at What Are the Types of Deconvolution, Sharpen Defocused Image (Deconvolution / Image Restoration), Deconvolution by Convolution and 1D Deconvolution with Gaussian Kernel (MATLAB).
Salt and Pepper Noise
Salt and Pepper noise is usually treated with Median Filter.
In the case of Color what you can do is either work on each channel on its own or generalize the cost function the Median minimizes to vectors (Namely minimize the $ {L1} $ Norm of the neighborhood of RGB vectors.
Contrast
The last degradation seems to be some kind of contrast removal from the image.
In order to reverse is normalize the image into the range $ [0, 1] $ (Black to White). You can do that by the following method:

Find the Minimum Value of the image, lets' call it minVal.
Find the Maximum Value of the image, let's call it maxVal.
Set the output value of the $ i, j$ pixel as following: mOuptutImage(i, j) = (mInputImage(i, j) - minVal) / (maxVal - minVal).


Sampling H(z) to get DFT Suppose that I have a $H(z)$ and I sample it to get a DFT of 15 values. Let's call this DFT $H_{1}[k]$. Then, suppose I antitransform $H(z)$ and grab the first 10 values of the sequence, and then I add 5 zeros. After this, I calculate the DFT of this new sequence. Let's call it $H_{2}[k]$. What differences would these two DFTs have between each other (if any)?

When grabbing the first 10 samples of $h_1[n]$ you are essentially multiplying the whole sequence $h_1[n]$ with the rectangular window

$$
  w[n] =
\begin{cases}
  1 & 0 \leq n < 10 \\
  0 & \mbox{otherwise}
\end{cases}
$$

whose DFT is given by:

$$
\begin{align}
  W\left[k\right]
    &= \sum_{n=0}^{9} e^{-j 2\pi k n / N} \\
    &= \frac{1 - e^{-j 2\pi k \cdot 10/15}}{1 - e^{-j 2\pi k / 15}}
\end{align}
$$

By the circular convolution theorem:

  Multiplication in the discrete-time domain becomes circular convolution in the discrete-frequency domain. Circular convolution in the discrete-time domain becomes multiplication in the discrete-frequency domain.

In other words,

$$
  \mbox{DFT}\left\{h_1[n] \cdot w[n]\right\}
  = \frac{1}{N} H_1\left[k\right] \otimes W\left[k\right]
$$

where $\otimes$ is the circular convolution operator.
This can be seen with the following matlab script (you can play around with h1 to see what happens with other waveforms):

n = [0:14]; % discrete-time index
k = [0:14];
h1 = sin(pi*n/14); % some arbitrary function

% window function
w = [ones(1,10) zeros(1,5)];

% time domain operation
h2  = h1.*w;
% convert to frequency domain
H2t = fft(h2);

% frequency domain operation
W   = fft(w);
H1  = fft(h1);
H2f = cconv(W,H1) / 15;

% Compare H2t & H2f
figure(1);
subplot(2,1,1);
hold off; plot(k, abs(H2t), 'b');
hold on;  plot(k, abs(H2t), 'bx');
hold on;  plot(k, abs(H2f), 'ro');
title('Magnitude');
subplot(2,1,2);
hold off; plot(k, angle(H2t), 'b');
hold on;  plot(k, angle(H2t), 'bx');
hold on;  plot(k, angle(H2f), 'ro');
title('Phase');

Note that if $h_1[n]$ is already no longer than 10 samples (i.e. if those truncated samples are already zero), then windowing in the time-domain will have no effect and similarly the circular convolution with $W[k]$ in the frequency domain will happen to produce $H_1[k] = \frac{1}{N} H_1[k] \otimes W[k]$ (but only as a special case).


reconstruction filter - How does it actually work? I'm trying to form my own understanding on the religious war around using 192kHz as a sampling rate for playback (the Internet seems to have a wealth of material on both sides). I'm struggling to understand how reconstruction filters work.

The Nyquist-Shannon sampling theorem, often cited by the anti-192 camp, basically states that a 44.1kHz sample rate is sufficient to reconstruct without loss a band-limited signal of 20kHz. However, when looking at the Whittaker–Shannon interpolation formula, it seems to me that an ideal reconstruction filter would need to have access to all samples, i.e., all past and all future samples.

I'm not an expert in analog audio, but I doubt that such a device can be built. I would assume that, at best, a delay can be introduced to "wait" for enough future samples to arrive, so that the contribution of unavailable future samples to the current output instant becomes negligible.

Can someone explain how practical reconstruction filters work and what are their tradeoffs? Is there a theoretically stricter limit to the Nyquist-Shannon theorem if only a window of samples is available or if reconstruction latency is unacceptable?

The sampling theorem requires a perfectly bandlimited signal, bandlimited to below twice the sampling frequency.  The problem with this is that only an infinite length signal (e.g. exists before the big bang) can be perfectly bandlimited.  This is from the Fourier theorem regarding any domain with finite support.

Thus all real-world signals are imperfectly band-limited, limited by the length of the recording session and the duration of the impulse response of the low-pass filtering used before sampling (e.g. imperfect band-limiting).  So you must assume a finite noise floor due to sample aliasing.  Therefore your reconstruction also need be no better than this noise floor.  Therefore you can window the perfect reconstruction formula to some reasonable finite duration.  For audio, this window width can be limited to on the order of a small multiple of the period of the lowest frequency audible, as stuff before and after that can't affect audible pitch perception significantly (other stuff, such as masking and adaptive thresholding, etc., likely dominate).

The advantage of very high sample rates is that the aliased frequencies (due to imperfect pre-sampling filtering and post-sampling reconstruction) are even higher, and far less likely to be picked up by a physical microphone or pass with any significant magnitude (above the thermal noise floor, etc.) through the anti-aliasing filters.  As other answers here detail, this also allows physically realizable filters to be flatter and have more linear phase response in generic 20-20k pass-band, and roll-off more approaching half the sample rate.


Intuition about tag pose estimation accuracy I have some difficulties to have a correct intuition about expected pose accuracy when dealing with planar pose estimation (e.g. pose estimation from AprilTag detection) and perspective camera model.
Assuming ideal conditions (e.g. simulated data), I consider also that the tag corners can be extracted with a good accuracy, e.g. < 1px.

Given two physical tags projected onto the image plane with the same size, oblique view (tag parallel to the image plane), can we assume the expected pose error is the same regardless of the working distance?

Since for tag detection, what matters is the size of the tag in the image, and for pose estimation what matters is precise corners location extraction, I want to put a minimal tag size in the image to choose other factors like physical tag size, camera focal length, working distance, ...
If 1. is true, this means that for two tags projected with a 100 px size in the image, if we have a translation pose error of 1mm at 1m, we could assume to have 1mm translation pose error at 2m if we choose appropriately the physical tag size to ensure a projected tag size of 100px.
My intuition for 1. is that it is not true, and somehow the translation pose error is dependent to the working distance.

For a fixed working distance, if a tag is projected onto the image plane with a 100px size, what is the impact of the projected tag size wrt to the translation pose error?

For instance, if the size of a tag is 100px in the image and the translation pose error is 1mm, if we double the size of the tag in the image to be 200px by playing with the focal length or by increasing the physical tag size, can we expect to reduce the translation pose error by the same factor (i.e. 0.5mm)?
For sure bigger tag size should improve the accuracy of the pose estimation, but I am not sure it scales like that.

By pose estimation, I mean the translation ($t_x, t_y, t_z$) + orientation ($3 \times 3$ rotation matrix) for rigid body transformation.
Also called PnP pose estimation, like pose from homography for planar object or EPnP by Lepetit et al.
For perspective projection:

normalized image coordinates
$$
x = \frac{X_c}{Z_c} \\
y = \frac{Y_c}{Z_c}
$$
to the image coordinates
$$
u = x \times f_x + c_x \\
v = y \times f_y + c_y
$$

For focal length in pixel, I have:
$$
f_xpx = \frac{f_xmm \times \text{image_width}}{\text{sensor_width_mm}}
$$

This intuition would be very useful to accordingly choose the physical tag size, the focal length, the image resolution to ensure or at least maximize the probability to have a given translation pose error.

The size of the tag is an important parameter. But it is a secondary parameter. The key parameter that we are looking at here (and one that is not easy to estimate) is Signal to Noise Ratio (SNR). So, very very briefly: The bigger the tag, the more data we obtain from the tag, the better the position of the estimate will be.
But, the "size of the tag" as perceived through a camera is an entirely relative concept. If you take a picture of a lighthouse that is far from the coast with a web cam at 320x240, you will see (for example) a flickering pixel patch of size $6 \times 2$. You can attach a lens to the camera or crank its resolution up. Both of these solutions, will allow you to image the lighthouse much better. If either of these options are not available then, there is a maximum distance (which puts a limit on the perceived size of the lighthouse) within which the lighthouse is sufficiently sampled.
So, the perceived size of the tag is only relevant with respect to the SNR. In low noise, small patches will be extracted more accurately and vice versa.
A discussion about accuracy must go through the way the detector derives the position of the tag.
Very briefly: There are two parts to it, detection (What pixels make up a "tag" in an image?) and decoding (Which tag(s) appear(s) in an image?). We leave decoding aside for a minute. The first step in detecting where the tag is is the computation of the gradient (magnitude and direction). The second step is the grouping of similar value gradients. Gradients (slopes) of similar value help in picking out lines. So each one of these groups (potentially) contains a line. We take the $(x,y)$ data of those groups and we use least squares to fit lines. At that point detection concludes. If you do this for all connected components, you will get back a data structure of connected lines that describes where the tag is (within the camera's field of view) and its orientation. There are some hints that we take from detection that are then used to decode the payload of the tag but this highlights why the tags are composed of black and white patches.
The simplest extreme counter-example that shows that the size of the tag is not really the issue here is this: Put a large enough tag in front of the camera and turn the contrast all the way down. This "kills" the gradient which is the main way by which the vision algorithm makes sense where the tag is. If the image is a grey fog, no tag will be detected because the gradient will return zero groups of pixels.
Conversely, suppose that the tag is of sufficient size, the contrast is optimal, the lighting is perfect, but, the code is amongst a huge herd of zebras chilling out in a field. That is, the background is now overwhelming the sensor with a huge amount of possible gradients that form lines that COULD lead to a tag just by chance.
So, the size of the tag is relevant with respect to the disturbance present in the image.
The noise, spreads out the gradient magnitudes and directions, which introduces uncertainty to the least squares of line fitting, which introduces uncertainties to the determination of position and orientation of the tag. The tag is black and white because this maximises the response of its detector. The gradient. It improves its chances of picking it out in the image.
We can derive some simple limits here of course. For example, you need at least 2 points per line segment of a tag so that all line segments that compose the tag have a chance of making it through the detector succesfully and the whole tag gets a chance to be recognised. So, the smallest perceived size for the tag is its size as if it was a set of pixel patches times 2. If you have a 4x6 tag (that is 4x6 black and white pixels), it needs to be imaged at least at 8x12 pixels. If not, you don't have 2 points per line to define it.
Then, what is the proportion of light that falls on the sensor that is not due to the tag? The variation of that light, blurs the gradient, which blurs the position of the recognised lines. Once the gradient is blurred the position is uncertain. The point is this, if you want to keep the ratio of signal created by the tag to signal NOT created by the tag constant, you need to vary the size of the tag so that you start enjoying the benefits of least squares that follows the gradient. More data, more chances of fitting the true line. Varying the size of the tag means either making it bigger, imaging at closer distance, imaging it through a larger magnitude lens, increasing the contrast, etc.
You might notice that this is not addressed in the original paper too. The authors provide results for recognition error rate versus distance and view angle but these are based on simulated data and also, they are not absolute (they cannot be absolute). They are comparative. They make sense with comparison to another method. (And also, they include the decoding of the tag content which puts additional resilience to the SNR we talked about before. While we are at this, you might want to see section "VI.b. Localisation Accuracy" and maybe section "III.Detector" which describe how the code is picked up in more detail than I go into here)
So, having said this, let's revisit the question:

"Assuming ideal conditions (e.g. simulated data), I consider also that the tag corners can be extracted with a good accuracy, e.g. < 1px."

If the conditions are ideal, the gradient-->least squares works perfectly and the pose is recognised to the limit of numerical accuracy (and image resolution).

Given two physical tags projected onto the image plane with the same size, oblique view (tag parallel to the image plane), can we assume the expected pose error is the same regardless of the working distance?

Yes, provided that the tags are within the limits of detection we talked about before. Side note: The tag that provides the most information to the gradient-->least squares part of the detector will have the best accuracy (between the two we are looking at).

For a fixed working distance, if a tag is projected onto the image plane with a 100px size, what is the impact of the projected tag size wrt to the translation pose error?

The question is ill posed. The size of the tag is only half the story. What else is going on around the tag?
I think that it will be very difficult to arrive at a quantitative estimation of this in a way that it would drive the selection of the tag size without some sort of experimentation or some sort of regularisation or control for some of these factors.
I hope this helps.


Do we have to insert the DC null in OFDM when using RF direct sampling ADC architecture? Generally, DC Nulls are added in OFDM because zero-IF receivers don't pass the DC component. But what if I use an RF direct sampling architecture? Will there be DC offset? If not, then is it okay not to insert a DC null?
Thank you in advance!

RF direct sampling architecture?
Will there be DC offset?

No, since the center of your signal doesn't end up on DC, so there's no DC offset.

If not, then is it okay not to insert a DC null?

yes.
However, you might want to use quadrature mixers for the other end, so you might not want to do that.
Also note that instead of direct sampling, superhet / low-IF approaches can also get rid of the DC offset, as does, to a not-so-bad degree, calibration. The question whether you need to insert a zero carrier at the center frequency of your OFDM signal is often a complex, more systemic one than just looking at the receiver architecture.


Generalized Barker sequences So we're doing a project where we've decided on QAM-n modulation. I've read up a little bit on synchronization and found Barker sequences, which supposedly has nice autocorrelation properties which makes them good for synchronization.

So one option is to use BPSK modulation only in the synchronization phase and then move onto QAM for the data. But to me, it's very unintuitive that we couldn't construct even NICER synchronization sequences (wrt autocorrelation peakedness etc) if we utilize more of the signal space. Maybe something like simultaneous Barker BPSK for both I and Q?

Is there a nice explanation for why it's generally best to use BPSK for the synchronization?

Longer preambles make it easier for the demodulator to lock onto the signal well.  This is due to getting more preamble signal energy and because longer preambles make it easier to detect stuff like carrier offset.

Barker sequences are nice for all the reasons you state, but they are also very limited.  The longest Barker sequence is 13 bits and there is only one of them.  The longest Barker sequence length that has two sequences is four, so you'd be limiting yourself to a four bit long preamble, which is not a very useful preamble.

You may say to yourself, "I could just use the same 13-bit Barker sequence for both the I and Q."  If you do that, though, you have made a BPSK constellation, only you've rotated it by $\frac{\pi}{4}$.

If you want to use a QPSK preamble, then I would use pseudo-random sequences from linear feedback shift registers (LFSR's).  You can find the polynomials for maximal length LFSR's here.