Multiplying signals in discrete-time vs continuous-time

You reached a puzzling conclusion about $c_1(t) = c_2(t)$, and wonder whether you made a mistake in deriving them, or if the equality is indeed correct then how to explain it, perhaps by explicitly deriving one from the other.
I cannot tell whether it's possible to explicitly manipulate the double summation in $c_2(t)$ so as to convert it into the single summation of $c_1(t)$. It may be possible, I haven't tried. But sometimes it's not possible to explicitly show it, and instead, we have to rely on indirect evidence to confirm the equality and use the equality for our advantage if possible.
One such example is the Fourier transform of the constant $1$ :
$$ \mathcal{F}\{1\} = 2\pi \delta(\omega) \tag{0} $$
The equality in Eq.0 is not derived by explicitly evaluating the forward Fourier integral, but instead, deduced from  the duality property of CTFT, and given that Fourier transform of $\delta(t)$ is $1$.
At the end of the analysis, we conclude that the validity of the equality $c_1(t) = c_2(t)$ is a consequence of the Nyquist sampling theorem;i.e., the truth of the equality is imposed by the sampling theorem, rather than a result of explicit algebraic manipulations of $c_2(t)$ into $c_1(t)$ or vice versa. And indeed, this is a useful side application of the theorem to prove that some equation is true which is very hard, if not impossible, to do so otherwise.
Let me show you, therefore, an indirect way of imposing the equality.
Let all signals $a(t),b(t)$, and $c(t)=a(t)b(t)$ are sufficiently bandlimited so that we can avoid aliasing.
Observe impulse train modulation relation:
$$ x_s(t) = x(t) \cdot \delta_T(t) ~ \cdot \cdot \cdot ~ \delta_T(t) = x(t) \cdot \delta_T(t) \tag{1}$$
where  $ ~\delta_T(t) = \sum_n \delta(t-nT)$.
Also observe the interpolation relation :
$$ \left( x(t) \cdot \delta_T(t) \right) \star h(t) = x(t) \tag {2}$$
where $~h(t) = \text{sinc}(t/T) ~$ is the ideal lowpass brickwall interpolation filter.
We will use Eqs. 1 & 2 to derive alternative but equivalent expressions for the samples $c_s(t)$ of $c(t)$ and achieve the single and double sum versions $c_1(t)$ and $c_2(t)$.
The Double Sum :
$$
\begin{align}
c(t) &= c_s(t) \star h(t) \\
&= \left( c(t) \cdot \delta_T(t) \right) \star h(t) \\
&= \left( a(t) \cdot b(t) \cdot \delta_T(t) \right) \star h(t) &\text{Multiply then sample}\\
&= a(t)b(t)  &\text{by Eq.2} \\
&= \left( \sum_n a[n] h(t-nT) \right) \left( \sum_m b[m] h(t-mT) \right) \\
&= \sum_n \sum_m a[n] b[n] ~ h(t-nT) ~ h(t-mT) \\
\end{align}
$$
The Single Sum :
$$
\begin{align}
c(t) &= c_s(t) \star h(t) \\
&= \left( a(t) \cdot b(t) \cdot \delta_T(t) \right) \star h(t) \\
&= \left( a(t) \cdot \delta_T(t) \cdot b(t) \cdot \delta(t) \right) \star h(t) &\text{by Eq.1} \\
&= \left( a_s(t) \cdot b_s(t) \right) \star h(t) &\text{Sample then multiply} \\
&= \left( \sum_n a[n] \delta(t-nT) \right) \left( \sum_m b[m] \delta(t-mT) \right) \star h(t) \\
&= \left( \sum_n \sum_m a[n]b[m] \delta((m-n)T) \cdot \delta(t-mT) \right) \star h(t) \\
&= \sum_n a[n] \left( \sum_m b[m] \left[ \delta((m-n)T)\delta(t-mT) \star h(t) \right] \right) \\
&= \sum_n a[n] \left( \sum_m b[m] \delta((m-n)T)) h(t-mT) \right) \\
&= \sum_n a[n]b[n] h(t-nT)
\end{align}
$$
The sampling relations are  :
$$a_s(t) = a(t) \delta_T(t) = \sum_n a[n] \delta(t-nT) \tag{3}$$
$$b_s(t) = b(t) \delta_T(t) = \sum_m b[m] \delta(t-mT) \tag{4}$$
$$c_s(t) = c(t) \delta_T(t) = \sum_n c[n] \delta(t-nT) \tag{5}$$
We can get back $a(t),b(t)$,and $c(t)$ by ideal bandlimited interpolation of their samples $a[n],b[n]$,and $c[n]$ :
$$a(t) = a_s(t) \star h(t)= \sum_n a[n] \text{sinc}((t-nT)/T) \tag{6}$$
$$b(t) = b_s(t) \star h(t)= \sum_m b[m] \text{sinc}((t-mT)/T) \tag{7}$$
$$c(t) = c_s(t) \star h(t)= \sum_n c[n] \text{sinc}((t-nT)/T) \tag{8}$$


Obtaining power spectral density without having knowledge about frequency bandwidth or sampling frequency of the signal in Matlab?

Without more information about the signal, it is not possible.

The set of samples is just a set of numbers with no further information associated with them. If you know the frequency at which the samples were taken, then you can associate "real-world" frequencies to their FFT bins. Otherwise, it may be impossible.

To give an example, consider the samples of $\cos(2\pi 100t)$ taken at $f_s=500\,\text{Hz}$, and the samples of $\cos(2\pi 200t)$ taken at $f_s=1000\,\text{Hz}$. They are exactly the same! There is no way to tell the two signals apart just from their samples.

It may be, however, that you have some extra piece of information that may help deduce $f_s$. For example, you may know that the signal's spectrum has a peak, or a particular phase, or some other property at a certain frequency. If you find the normalized FFT of the samples and identify this property, you can extrapolate from there to estimate the actual signal frequencies.


What happens to the audio of a reconstructed signal if we remove too much frequency?

Theoretically, the error between the original and the reconstructed signal goes up.  Practically, for actual audio, if you're lowpass filtering too much, the audio starts to sound muffled or muddy -- like you're trying to listen to the event through a pillow.
Again in theory, if you've got a bandlimited signal mixed with noise, the best processing for that signal is to low-pass filter it more than the signal is already bandlimited.  There's a lot of theory out there on how to design optimal filters for this.


Can we generate signals in reality that have non symmetric magnitude spectrum in the Fourier domain?

As @Bob said, it depends of your definition of "in reality". I am going to ignore images on this post and focus on signals as a function of time.

One-dimensional signals always have a symmetrical magnitude spectrum. Examples are the voltage or current on a wire, or a sequence of real numbers stored in digital memory.

Two-dimensional signals may have a non-symmetrical magnitude spectrum. Examples are complex signals, whether defined mathematically or stored in a digital memory; and the voltage or current on two wires, one of which is labeled $\sqrt{-1}$ and with the agreement of doing complex arithmetic (see $9\rm V$ Battery with $45^\circ$ phase).

In my opinon, all these signals exist in reality, and then the answer is that yes, some signals in reality may indeed have a non-symmetrical magnitude spectrum.


Overall system IR: How to do arithmetic with $\delta$ functions?

You make it unnecessarily complicated.

Given
$$\begin{align}
h_1[n]&=2\delta[n-2]-3\delta[n+1]\\
h_2[n]&=\delta[n-1]+2 \delta[n+2]
\end{align}$$
to calculate
$$h[n]=h_1[n]\star h_2[n]$$

the easiest way is to use the sampling property of delta, which is

  $$f[n]\star\delta(n\pm n_0)=f[n\pm n_0]$$

So you can write as follows

$$\begin{align}
h_1[n]\star h_2[n] &= h_1[n]\star \left(\delta[n-1]+2 \delta[n+2]\right)\\
&=\color{green}{h_1[n-1]}+2\color{blue}{h_1[n+2]}\\
&=\color{green}{2\delta[(n-1)-2]-3\delta[(n-1)+1]}+2\color{blue}{\left(2\delta[(n+2)-2]-3\delta[(n+2)+1]\right)}\\
&=2\delta[n-3]-3\delta[n]+4\delta[n]-6\delta[n+3]
\end{align}$$


As far as I know there are no known tight upper and lower bounds, but there's a quite accurate empirical formula for estimating the required filter order $M$ given maximum passband and stopband deviations $\delta_1$ and $\delta_2$, and the width of the transition band $\Delta\omega$:
$$M=\frac{-10\log_{10}(\delta_1\delta_2)-13}{2.324\;\Delta\omega}\tag{1}$$
Formula $(1)$ can be found as Eq. $(7.117)$ in Oppenheim and Schafer's textbook Discrete-Time Signal Processing (3rd ed, p.570). According to the authors, this formula is due to Kaiser, who developed it based on experimental results published by Herrmann, Rabiner, and Chan.
Also Matlab's firpmord.m uses a similar but more complex formula to estimate the required filter order.


The OP stated he was interested in $\pi/4$-DQPSK (not QPSK), so phase synchronization is presumably not an issue for him.

As far as the actual implementation is concerned, you'll save yourself some time if you become familiar with the bottom of pages 29 (symbol mapping) and 37 (differential detection) in this student paper. Ignore all the old TI DSP-chip specific information and the information on pulse-shaping, which is already built into GNU Radio (if the OP is using that). There is also a reasonable description of a simple timing recovery algorithm. It is serviceable in some scenarios, but not others.


The optimality criterion for Bessel (low pass) filters is that they have a maximally flat group delay at $\omega=0$, as correctly pointed out in a comment by Peter K.. The derivation is slightly tortuous, but you can find it for the continuous-time case in these lecture notes, which are taken from the book Introduction to Circuit Synthesis and Design by Temes and La Patra. I can recommend this book to anybody who wants to gain a deep understanding of electrical network theory.

The interesting fact is that unlike for other standard filters (Butterworth, Chebyshev, Cauer), the discrete-time version of a Bessel filter cannot be derived from its continuous-time counterpart via the bilinear transform. The reason is the frequency warping inherent in the bilinear transform, which doesn't affect the optimality conditions of the other standard filters (maximal flatness of the magnitude, or equiripple behavior of the magnitude, or a combination thereof), but it does not preserve the maximally flat group delay. A discrete-time filter derived from an analog Bessel filter via the bilinear transfer will not be optimal in any way, and it will not have a maximally flat group delay. This means that a discrete-time Bessel filter must be derived directly in the discrete-time domain. This has been done in this famous paper by Thiran. Again, the derivation is slightly painful but I'm afraid there is no easier route.

What you can take away:

Bessel filters have a maximally flat group delay.
Discrete-time Bessel filters cannot be obtained from continuous-time Bessel filters via the bilinear (or any other) transform.


RF mask (power spectrum measurement bandwidth)

The simple answer: when you take a power measurement using some piece of test equipment (such as a spectrum analyzer or an RF power meter), the reported value indicates the amount of power measured within some detection bandwidth of the instrument. For some devices, that might be a fixed value, determined by the filtering present in its RF front end. Other types of equipment, like a spectrum analyzer, will typically have adjustable bandwidth (which in this case would usually be called the instrument's resolution bandwidth).

To see how the bandwidth setting affects the measurement, try injecting a modulated signal into a spectrum analyzer and see how the power measurements vary as you change the resolution bandwidth: you can get very different answers. Take a common $\frac{\sin{x}}{x}$ spectral envelope as an example:

If you use a very small measurement bandwidth, then you're taking power measurements over very small swaths of frequency. This allows you to observe the detail in the spectrum's shape.
As you increase the resolution bandwidth, your instrument will take more and more of the signal's spectrum into account for the power measurement, and the reported values will become more of an overall measure of the total received power.

As was indicated in the comments already, the specification calls out a specific bandwith to ensure repeatability. Specifications are used to validate implementations as conforming to an existing standard. This particular scheme for measuring the emission mask of a transmitter implementation is likely in use because that was the one that was chosen when the standard was designed. Due to the imprecise definitions of concepts like SNR and bandwidth, there are many ways that you could define such a requirement, but this is the one that the authors chose. That doesn't make it right or wrong; it just is.


Design a comb notch filter when sampling rate is not an integer multiplicative of the cutoff frequency

You can design that filter manually without problems. Matlab just uses a very simplistic approach to comb filtering with a delay line.

In order to keeps things as simple as possible I would recommend you use a series of notch filters to remove each partial of your harmonic noise separately. That also gives you more control over how much of each harmonic really has to be removed.

A simple notch filter is one with unit gain at DC, i.e. $H(1) = 1$, and the zeros on the unit circle at the desired notch frequency, positive and negative. The poles have to be inside the unit circle, close to the zeros to cancel their effect farther away from the notch frequency. We can design this directly in the z-domain.

$$H(z) = \frac{(z-\exp(i\omega))(z-\exp(-i\omega)}{(z-r\exp(i\omega))(z-r\exp(-i\omega))}\times\frac{(1-r\exp(i\omega))(1-r\exp(-i\omega)}{(1-\exp(i\omega))(1-\exp(-i\omega))}$$

Here $\omega$ is the normalized frequency, meaning $\omega=2 \pi \frac{f_0}{f_s}$ where $f_0$ is your notch frequency and $f_s$ the sampling frequency. The parameter $r\in]0,1[$ controls the width of the notch. The second fraction makes sure we get a unit response at $\omega=0$.

To get the coefficients for a recursive filter implementation, we have to simplify the transfer function and cancel $z^2$ to turn it into a rational function of $z^{-1}$. The result is
$$H(z) = \frac{\left(1+r^2-2r\cos(\omega)\right)\left(z^{-2}+1-2z^{-1}\cos(\omega)\right)\csc(\omega/2)^2}{4\left(r^2 z^{-2}+1-2rz^{-1}\cos(\omega) \right)}$$

So you can read off your filter coefficients to be

$$ A[n]=\left(1,-2r\cos\omega,r^2 \right)$$
$$ B[n]=\left(1,-2\cos\omega,1 \right)\cdot\frac{1}{4}\left(1+r^2-2r\cos\omega\right)\cdot\csc(\frac{\omega}{2})^2$$

The last open question is how exactly the parameter $r$ controls the notch behavior. Near the zero the transfer function on the unit circle locally behaves like

$$f(x) = \frac{x}{x-i(1-r)}$$

around $x=0$. For $r=1$ we get a singularity at $x=0$ and an ill-defined transfer function. For values $r<1$ but close to $1$ we can see that for both $x\to-\infty$ and $x\to\infty$ the transfer function becomes $1$. The interesting question is, how far does the effective influence of the pole reach on the x-axis?

The square magnitude of $f(x)$ is $$|f(x)|^2=f(x)\cdot f(x)^*=\frac{x^2}{x^2+(1-r)^2}$$
which is always real and positive. We can easily equate it to a threshold and solve for $x$.The most common threshold for filter cutoffs in signal processing is $1/\sqrt{2}$, and if we solve $|f(\Delta x/2)|^2==1/\sqrt{2}$ we get the effective bandwidth of the notch to be
$$\Delta x = 2(1-r)\sqrt{1+\sqrt{2}}$$

and we can use the equivalence of our problem to the original pole/zero placement to conclude

$$ r = 1-\frac{\Delta \omega}{2\sqrt{1+\sqrt{2}}} $$

with the natural frequency bandwidth $\Delta \omega$. If you want this as a true frequency then again $\Delta \omega=2\pi\frac{\Delta f}{f_s}$

Note that this bandwidth calculation used a local approximation, so it is strictly only true for $(1-r)<<1$, or small bandwidths compared to the sampling frequency.

So this gives you a notch filter for one single harmonic of your harmonic noise. Just apply one filter at a time, one after the other, with the parameters matched to the harmonic you would like to cancel, up to fs/2. The matlab filter() function readily takes the coefficients A[n] and B[n] given above.

An attempt to implement this in scilab seems to work nicely.  Plot of before and after spectrum:

Code:

// 12740
fs = 1000;
f = 101.89798;
omega = 2*%pi*f/fs;
delta_omega = omega/2;
N= 1000;
t = [0:N];
f = [0:N]/N*fs-fs/2

r = 1 - delta_omega/(2*sqrt(1+sqrt(2)));

num = [1 -2*cos(omega) 1]*1/4*(1+r.^2-2*r*cos(omega))*csc(omega/2).^2;
den = [ 1 -2*r*cos(omega) r.^2];

X = sin(omega/2*t) + sin(omega*t) + sin(2*omega*t);
Y = filter(num,den,X);

figure(1)
clf;
plot(f,log(fftshift(abs(fft(X)))+0.01))
plot(f,log(fftshift(abs(fft(Y)))+0.01),'k:')


What you describe is a direct-sampling architecture. Those are practically available for radio amateurs using frequencies up to 30MHz. They are also actively researched to open up their usage for much higher frequencies. Many interesting topics (e.g. using the aliasing effect to your advantage and undersample the signal) are investigated here and the possibilities of these receivers are large. However, it will probably take a long time until those are practically available for the WiFi frequency range.

However, your WiFi transceiver is most likely a direct-conversion architecture. In the digital domain of it, the data-carrying signal is generated. In case of 802.11n it is an OFDM modulated signal with up to 20/40MHz bandwidth. This signal is centered around 0Hz and hence cannot be transmitted as it is. This signal is called the baseband signal.

After Digital/Analog-conversion, in a stage called a mixer, this signal is then upconverted to the final carrier frequency around 2.4GHz. In the spectrum, the signal is literally shifted from a center frequency of 0Hz to a center frequency of 2.4GHz. At this point, the signal is called a passband (or bandpass) signal. After amplification, this signal is then transmitted via the antenna.

On the receiver side, the same happens the other way around. There is another mixer that downconverts the passband signal from 2.4GHz center frequency to a baseband signal at 0Hz center frequency. At this point, the signal has just a bandwidth of 20/40MHz and can easily be sampled in an A/D converter.

What follows is demodulation, error correction and other digital processing to eventually show you pictures of cute kittens in your browser.


Advice on designing a digital filter that doesn't have phase-sensitive edge artifacts?

Your basic problem is that filtfilt (and most other linear filtering routines) take filters that are designed for infinitely long time expanses, and apply them to a chunk of data as if the data were extended infinitely in both directions with zeros.
So you have a legitimate bandpass filter, and it's "seeing" a legitimate jump in the signal at the starting point of your signal.
There's three basic approaches you can take; the first two are ad-hoc and easy, the third one is difficult if you're starting from first principles.  It's certainly been solved out there someplace, but a brief search here on "filter finite-length data" didn't find me joy.
Approach 1: window the input data
Take your input data, and multiply it by something that'll make it taper off at the ends.  I.e. a ramp from 0 to 1 over 10 samples at each end, or $\frac{1}{2}\left (1 - \cos \frac{\pi n}{N} \right)$ for N samples at each end (suitably reversed on the left end).  You'll have some artifacts (a rising sine wave isn't the same as a steady one, after all), but they'll be attenuated.
Here's python code implementing a cosine edge attenuation with the ability to customize what central % of the signal is kept as 1:
def attenuate_edges(signal,time,edge_attenuation_percent):
  start = int(np.floor(len(time)*edge_attenuation_percent))
  end = int(len(time)-start)
  ramp = (1-np.cos(np.pi*(np.arange(start)/start)))/2
  edge_attenuator = np.ones(len(time))
  edge_attenuator[0:start] = ramp
  edge_attenuator[end:len(time)] = np.flip(ramp)
  return(signal*edge_attenuator)

Approach 2: Trim the output data
Do what you're doing now, and lop off the nastiness at the ends.  This is probably the easiest, and if you can just collect a bit more data, doesn't lose you anything.
Approach 3: Do a proper estimate of the interfering signal, and subtract it out
This will be fun-fun if you love math and have the time.  Basically you'll use the fact that the value of your interfering signal at time $n$ correlates in a specific way with the values of your interfering signal at time $k$ for all values of $n$ and $k$ in your data set.  You'll probably end up with something that looks a lot like a Wiener or a Kalman filter, that takes the end-effects into account.  Your estimate will be worse at the ends, but this will show up as a bit of noise on the ends -- not as those honkin' big pulses.
If I couldn't figure out the search terms for this, it would take me a day to do and another day to verify, and supposedly an expert.  OTOH, Gauss or Laplace probably invented it in the 19th century, and may even have thought it important enough to write down, somewhere.  So I'm sure the method exists.


I have a simple compensation approach that I've used to implement a reasonable estimate of an inverse $\textrm{sinc}$ for use as a CIC compensator as well as other inverse $\textrm{sinc}$ applications. This approach makes use that the $\textrm{sinc}$ function in the passband can be reasonably approximated by a weighted cosine function. Therefore a raised cosine function in frequency, meaning $b-(\alpha)cos(\omega)$ can be used in cascade with the CIC filter as compensation. This function is adjusted by changing the cosine weight $\alpha$ while keeping $b=1+\alpha$, such that the mean squared error over the passband of interest of the cascaded result is minimized.

Once the minimized passband is established, (you can do this with least squared minimization techniques, but I have done it quickly by evaluating passband error for a least squared miminum while doing a binary search on $\alpha$ values, which quickly converges to an $\alpha$ value that minimizes the error for a given passband), the coefficients of the filter, which is the filter impulse response, are as follows:

Coeff 1, 3: $-\alpha/2$

Coeff 2: $b= 1 + \alpha$

It is that easy!  The filter is a simple 3 tap FIR with coefficients [$(-\alpha/2)$   $(1+\alpha)$   $(-\alpha/2)$]. Since the filter is symmetric, it can also be implemented with just two multipliers as shown in the figure. Symmetric filters also have the nice benefit that they are linear phase.

See the figure below for an example of using this CIC compensator in an x8 CIC Intepolator.  In use, the signal is passed through the 3 tap compensator at the lower rate, and then fed into the x8 CIC Interpolator with an output at 8x the sample rate.  The cascaded response shows the significant improvement that can be achieved with this simple compensator.

What should now be clearer from this figure is that this compensation approach will slightly decrease overall rejection in portions of the stop band. (For this specific example the decrease was on the order of 1.6 dB). However, since the unity gain positions of the compensator will be aligned with the null locations in the CIC response, this reduction in rejection is not in those rejection regions that matter most for CIC rate conversion.

Here is a zoom in of the passband showing the excellent flatness that can be achieved with three taps!

Note on solving for the FIR coefficients from the first figure above.

The coefficients being the impulse response of the filter are determined from the inverse DFT of the frequency response. They can be quickly derived from the generic equation for the frequency response for an FIR filter in terms of the FIR coefficients as:

$H(\omega)= \sum_{n=0}^{N-1}c_ne^{-jn\omega}$

where $c_n$ are each of the $N$ coefficients, and $\omega$ is the digital frequency domain from $0$ to $2\pi$ corresponding to 0 Hz to the sampling rate.

As for the case of our 3 tap symmetric FIR this is

$H(\omega)= -\frac{\alpha}{2}+ be^{-j\omega}-\frac{\alpha}{2}e^{-j2\omega}$

$ = e^{-j\omega}(-\frac{\alpha}{2}e^{+j\omega}+b-\frac{\alpha}{2}e^{-j\omega})$

$=e^{-j\omega}(b-\alpha(\frac{e^{+j\omega}+e^{-j\omega}}{2}))$

$=e^{-j\omega}(b-\alpha cos(\omega))$

which is exactly the desired frequency response as given in the first plot (the first term $e^{-j\omega}$ is just the required delay for the filter to be causal but does not affect the magnitude response (and shows how the phase is linear versus $\omega$ ...the same reason all symmetric FIR filters are linear phase- the response for any symmetric FIR can all be described in terms of cosines with a linear phase delay pulled out of the equation!).

It may be more intuitive to some imagine this operation when going from the time domain to the frequency domain if you are familiar with the FT of a cosine wave: A raised cosine wave in the time domain (meaning a cosine with a DC offset) has a Fourier tranform with 3 impulse components in the frequency domain; a DC term, and a positive and negative frequency term. Likewise, the same waveform in the frequency domain will have three impulse components in the time domain, and these components are the impulse reponse by definition, which are therefore the filter coefficients.


How can the Poles of the Root Locus be negative?

If the denominator polynomial is given by

$$Q(s)=(s-s_0)(s-s_1)\ldots (s-s_{n-1})\tag{1}$$

then the corresponding time domain function is made up of a sum of scaled exponentials of the form $e^{s_it}u(t)$, where $u(t)$ is the unit step function. These terms decay only if $\textrm{Re}\{{s_i}\}<0$, so for stability the zeros of the denominator polynomial (i.e., the poles of the system function) must lie in the left half-plane.

Let $h(t)$ be the impulse response of a causal LTI system. The transfer function is the Laplace transform of $h(t)$:

$$H(s)=\int_{0}^{\infty}h(t)e^{-st}dt\tag{2}$$

This integral only converges in a region $\textrm{Re}\{s\}>c$ ("region of convergence"), with some real-valued constant $c$. That constant depends on the properties of $h(t)$. If $h(t)$ is a weighted sum of exponentials $e^{s_it}u(t)$ then $c$ equals the most positive real part of the complex numbers $s_i$:

$$c=\max_i\big\{\textrm{Re}\{s_i\}\big\}\tag{3}$$

This makes sure that the damping by the term $e^{-st}$ in $(2)$ is enough to make $h(t)e^{-st}$ decay sufficiently such that the integral converges.

For stability we require $c<0$, i.e., all poles lie in the left half-plane, and the imaginary axis of the $s$-plane $s=j\omega$ lies inside the region of convergence.


Why can't we have IoT from 2005? but begin now,or future?

IoT requires inexpensive low-power communication. There were no solution with those 3 characteristics in 2005. Also, it was harder to monetize data so the perception was that IoT would never payback.


How to generate IQ component for transmitting 2FSK

Your I & Q plots look wrong to me.  Instead of trying to sort out what you may have done wrong, let me just give you a working example instead.

The following GNU Octave (Matlab clone) script shows how to generate the I & Q of 2-FSK with square pulses and a modulation index of $1/2$ (which should correspond to a bit rate of $4\cdot\Delta f$).

The "magic" happens on the lines following the comment "Frequency modulate the baseband packet"

% Keep Octave happy
pkg load signal;

% Bit rate
Rb = 9600;

% Modulation index; 1/2 is Minimum Shift Keying.
modulation_index = 1/2;

% Optional frequency shift for demonstration of modulation on a
% carrier within -Fs/2 to +Fs/2
freq_shift_hz = 0;

% Samples per symbol and sample rate out of the modulation process
sps = 10;
Fs = Rb * sps;

% pulse filter taps - square pulses
taps_pf = ones(sps, 1)/sps;

% Create a random bit string
nbits = 16;
message_bits = round(rand(nbits,1));
nrz_message_bits = (message_bits - 0.5)*2;

% Upsample and pulse shape the packet
% N.B. Octave's filter function sets up a history of 0's for us.
x = sps*[upsample(nrz_message_bits, sps); zeros(1, 1)];
packet_baseband = filter(taps_pf, [1], x);

% Throw in a frequency shift (as a level shift of the unscaled baseband)
% for demonstration purposes
freq_shift = freq_shift_hz / (Rb/2 * modulation_index);
packet_baseband = packet_baseband + freq_shift;

% Frequency modulate the baseband packet
fm_gain = pi/(Fs/2) * Rb/2 * modulation_index;
x = packet_baseband * fm_gain;
phase = cumsum(x); % phase is integral of frequency
packet_modulated = exp(1i*mod(phase, 2*pi));

figure(1);
stem(message_bits);
title('Message Bits');
xlabel('Bit number');
ylabel('Value');
grid on;

figure(2);
N1 = length(packet_baseband);
t1 = [0:N1-1]/Fs;
plot(t1, packet_baseband * Rb/2 * modulation_index, 'x-');
title('Upsampled, Pulse Shaped, Scaled, and Level Shifted NRZ Bits');
xlabel('Time (seconds)');
ylabel('Frequency Deviation (Hz)');
grid on;

figure(3);
N2 = length(phase);
t2 = [0:N2-1]/Fs;
plot(t2, phase, '.-');
title('Accumulated Phase');
xlabel('Time (seconds)');
ylabel('Accumulated Phase (radians)');
grid on;

figure(4);
N = length(packet_modulated);
t = [0:N-1]/Fs;
plot(t, real(packet_modulated), t, imag(packet_modulated));
title('I & Q of 2-FSK Modulated Signal');
xlabel('Time (seconds)');
ylabel('Amplitude');
grid on;


Do symmetric discrete signals have zero phase?

Assuming you are referring to the standard definitions of odd and even symmetry, i.e.
\begin{align}
\text{even: }&& x[-n] &= x[n] \\
\text{odd: }&&x[-n] &= -x[n]
\end{align}
then it still holds that the DFT of an even sequence is real (i.e. the phase is 0 or $\pi$) and the DFT of an odd sequence is imaginary (i.e. the phase is $\pi/2$ or $-\pi/2$).
This is independent of whether the the length of the sequence is even or odd as long as the sequence meets the definition of symmetry. This implies that for odd sequences $x[0]$ must be $0$ and for odd sequences of even length $x[N/2]$ must be zero too, where $N$ is the DFT length.
Both sequences that you show in your original question are not "symmetric" in the mathematical sense. They are only symmetric when properly time shifted and time shifting is clearly a non zero-phase operation.
The following Matlab script illustrates that this is true within the numerical accuracy of double precision math:
%% script to confirm zero phase DFT of even/odd symmetric sequences

N = 128; % FFT length

%% Even length sequences

% create a random sequence of length N/2+1
x = randn(N/2+1, 1);

% create even version
xEven = [x; x(end-1:-1:2)];

%  perform DFT and check zero phase
fx = fft(xEven);
zeroPhaseError = 10 * log10(sum(imag(fx).^2) / sum(abs(fx).^2));
fprintf('Even symmetry, Even length: zero phase error = %6.2fdB\n', zeroPhaseError);

%% Odd length sequences

% create odd symmetry version. This requires the samples at n=0 and n=N/2 to be
% zero
xOdd = [0; x(2:end-1); 0 ; -x(end-1:-1:2)];
%  check zero phase
fx = fft(xOdd);
zeroPhaseError = 10 * log10(sum(real(fx).^2) / sum(abs(fx).^2));
fprintf('Odd Symmetry, Even Length: imaginary phase error = %6.2fdB\n', zeroPhaseError);

%% now try an odd length sequence
% create a random sequence of length N/2
x = randn(N/2, 1);

% create even symmetry version
xEven = [x; x(end:-1:2)];

%  check zero phase
fx = fft(xEven);
zeroPhaseError = 10 * log10(sum(imag(fx).^2) / sum(abs(fx).^2));
fprintf('Even symetry, odd length: zero phase error = %6.2fdB\n', zeroPhaseError);

% create odd symmetry version. This requires the samples at n=0 to be zero.
% There is no sample at n=N/2 since that's not an integer.
xOdd = [0; x(2:end); -x(end:-1:2)];
%  check zero phase
fx = fft(xOdd);
zeroPhaseError = 10 * log10(sum(real(fx).^2) / sum(abs(fx).^2));
fprintf('Odd symmetry, odd length: imaginary phase error = %6.2fdB\n', zeroPhaseError);


changing filtering "speed" (cutoff frequency) of an IIR filter knowing only its coefficients

Your probably your best option here is "frequency warping". You can calculate the poles and zeros from the coefficients, warp the poles and zeros and then recalculate the coefficients again.
Warping is a procedure the applies a conformal mapping to the poles/zeros. Specifically a conformal mapping that maps the unit circle onto itself will maintain the overall shape of the filter but either bunch it up or stretch in the low frequencies (and vice versa in the high frequencies). A good choice is a first order allpass filter, that indeed maps the unit circle onto itself.
Loosely speaking you replace all delays in your filter with a first order all pass filter and recalculate the poles and zeros. That takes a bit of math to work out the details, but the actual code would be quite efficient.


Generate binary signal and convert to integer in Matlab

@Ali23. If by "How can I get a blue lines?" you mean how to generate samples of the your red curve, try the following:
m = [1, -1, -1, 1 1, -1]; m = [m m m m];
N = length(m);
M = 10; mup=zeros(1,NM); mup(1:M:NM)=m;  % oversample by M
h = hamming(M);                     % blip pulse of width M
x = filter(h,1,mup);     % convolve pulse shape with 'm' data
t = 1/M:1/M:length(x)/M;
figure(1), subplot(2,1,1), plot(0:0.1:0.9,h)
subplot(2,1,2), plot(t,x)
In the above I used a simple hamming sequence rather than the impulse response of a RRC filter. (You should use the impulse response of a RRC filter.) The key point here is that you have to upsample (zero stuffing) the 'm' sequence before you apply it to the filter.


Prewitt kernel is rotated 180 degree between correlation and convolution.

You can check this link: http://www.equalis.com/blogpost/731635/129540/Convolution-and-Correlation-in-Image-Processing


Notch filter: differences between IIR and FIR filters

Your assumption why IIR filters can have steeper transitions from passbands to stopbands compared to FIR filters of the same order is correct: IIR filter have poles away from the origin of the complex plane, and poles inside the unit circle close to zeros on the unit circle cause the corresponding frequency response to change rapidly with frequency.
The FIR filter in your question has a much broader notch than the IIR filter. This is not only caused by the fundamental difference between FIR and IIR filters, but also because for some reason the FIR filter has two zeros close to the notch frequency instead of only one. This makes the filter more robust with respect to errors in the estimation of the noise frequency, but it may also attenuate desired frequency components. It's also much harder for an FIR filter to approximate a constant response away from the notch frequency because the frequency response of an FIR filter is a polynomial and not a rational function as is the case for IIR filters.
From the zero locations of the FIR filter you can see that the filter is a linear phase filter: the zeros are either on the unit circle or they are mirrored at the unit circle. Note that FIR filters don't need to have linear phase. The required filter order for a certain specification can often be somewhat reduced if we don't impose the linear phase constraint.


Does $y[n] = x[n] \star (u[n]-u[n-2])$ have memory or is it memoryless?

The system $$y[n] = x[n] \star (u[n]-u[n-2])$$ where $u[n]$ is the unit step function, has memory.

Indeed the system is equivalent to $$y[n] = x[n] \star ( \delta[n] + \delta[n-1] ) \implies y[n] = x[n] + x[n-1]$$ and as it's clear from the given I/O relationship, the current value of the output $y[n]$, depends on the values input $x[n]$ at other times namely $x[n-1]$ (one step past value); hence the system has memory.


New to audio processing, why is my data so lopsided when I try to break into frequency ranges?

Most of the energy in an audio/speech signal is almost always found in the lower bands (roughly below 1kHz) so the lopsided shape that you're observing  is not surprising.

Let me just add some more information regarding fft frequency bins. These bins can be organized into larger buckets if you need this kind of representation.

For an FFT length of 1024, you should end up with 513 frequency bins (N/2+1). The first bin corresponds to DC (0 Hz) and is usually ignored and N/2+1 is the theoretical Nyquist frequency (also ignored).

The bandwidth of a frequency bin is defined as

BW_bin =  Sampling_rate/FFT length;

Note that while the sampling rate isn't necessary to compute the FFT, it is needed to calculate the bandwidth (frequency resolution).

To get the activity (magnitudes in dB), you can use the following equation (one way to compute the magnitudes)

Mag[i] = 10*log*(sqrt(2*(Real[i]^2+Img[i]^2/fftNorm));//iterate through bins [1-511]

fftNorm depends on the kind of the window function used (https://en.wikipedia.org/wiki/Window_function) and is simply N (fft length) in the case of a rectangular window (no window).
The factor 2 in the equation accounts for the upper (discarded) half of the FFT.

For visualization purposes, you can now easily combine several bins into a group of bins coressponding to a particular frequncy range.


Does spectral leakage include negative frequencies?

Yes spectral leakage includes both negative and positive frequencies; i.e, the convolution is from $-\pi$ to $\pi$.

However its effect depends on what exists at the negative frequency region, when you want to compute the resulting leakage on a given positive frequency.

On a typical application where you have a real valued sine wave of frequency $f_0$ as the windowed signal, then there is a dirac impulse at the negative frequency $-f_0$ as well as the positive frequency $f_0$. Then both impulses shift the spectrum of the window Fourier transform to their own locations and create leakage effects on any other frequency (unless they have zero value at that particular frequency). However the leakage is expected to be smaller from a negative frequency as it's side lobes are already decaying more than the positive frequency contributions. (This can be the opposite when both impulses are close to $\pi$ or close to $0$ radians.)

But if you have a complex signal, such as an analytic signal which has zero negative spectral content, then effectively the negative frequencies will not create leakage.


For a causal system to be stable, all the poles have to be inside the unit circle. In general, they may be complex, although they always come in complex conjugate pairs in order to keep the output real for real inputs.

So, the answer is no, there are no extra requirements for the poles to be on the real axis. However, they must be complex conjugates for the output to be real for a real input.

I would like to point out that the region of convergence of the Z Transform may include the unit circle, and yet the system may be non causal (but stable). For example, consider the impulse response: $h[n] = -a^{n}u[-n-1]$, for $|a| > 1$, which corresponds to a system having a pole outside the unit circle, and yet it is a stable (although non-causal) system.


Most certainly not. While there has been some claims to break Shannon here and there, it usually turned out that the Shannon theorem was just applied in the wrong way. I've yet to see any such claim to actually prove true.

There are some methods known that allow for transmission of multiple data streams at the same time on the same frequency. The MIMO principle employs spatial diversity to achieve that. Comparing a MIMO transmission in a scenario that offers high diversity with the Shannon limit for a SISO transmission in an otherwise similar scenario might actually imply that the MIMO transmission breaks Shannon. Yet, when you write down the Shannon limit correctly for the MIMO transmission, you again see that it still holds.

Another technique to transmit on the same frequency at the same time in the same area would be CDMA (Code Division Multiple Access). Here, the individual signals are multiplied with a set of orthogonal codes so that they can be (perfectly in the ideal case) separated again at the receiver. But multiplying the signal with the orthogonal code will also spread its bandwidth. In the end, each signal employs much more bandwidth than it needs and I've never seen an example where the sum of the rates was higher than Shannon for the whole bandwidth.

While you can never be sure that breaking Shannon is actually impossible, it is a very fundamental law that stood the test of time for a long time. Anyone claiming to break Shannon has most likely made a mistake. There needs to be overwhelming proof for such a claim to be accepted.

On the other hand, transmitting two signals on the same frequency at the same time in the same area is easily possible using the correct method. This is by no means an implication that Shannon is broken.


Frequency modulation implementation without trig functions

This is what an NCO (Numerically Controlled Oscillator) does. See this post for more details on the implementation and design considerations Numerically Controlled Oscillator (NCO) for phasor implementation?

By combining this with a D/A converter for an analog output we get what is commonly referred to as a Direct Digital Synthesizer (DDS).

With an NCO you can it you can easily do FM, PM and if you add an amplitude control at the output AM as well! To do FM specifically, you update the Frequency Control Word on each cycle according to the frequency vs time desired.

Since your question indicated a concern with efficiency; there has been a lot of work done in optimizing this approach to be highly efficient as well in terms of minimizing ROM storage. This includes more obvious approaches such as storing only 1/4 cycle of a sine wave, but also interpolation approaches are very common and other approaches include the Hutchinson Algorithm and Sunderland Algorithm.

The Hutchinson Algorithm uses a coarse ROM for a major phase step and fine ROM for a smaller increment and makes use of the small angle approximation:

$$\theta =\theta_F+\theta_C$$
$$sin(\theta) = sin(\theta_F)cos(\theta_C)+cos(\theta_F)sin(\theta_C)$$

for small anlge $\theta_F$:

$$sin(\theta) = \theta_Fcos(\theta_C)+ sin(\theta_C)$$

Thanks Hutchinson!

Also if you have time to spare in the implementation (meaning do not require a final update on each cycle), the CORDIC rotator is also quite interesting in terms of highly efficient approaches for implementing an NCO (or any other applications requiring phase rotation or trigonometric solutions).

References on the CORDIC:

Andraka, Ray, “A survey of CORDIC algorithms for FPGA based computers”, Proceedings of the 1998 ACM/SIGDA sixth international sympsosium on Field programmable gate arrays, pp191-200,  1998
Available at  http://www.adraka.com/files/crdcsrvy.pdf

The CORDIC FAQ at dspguru.com  http://dspguru.com/dsp/faqs/cordic

Maharatna, K, “Virtually scaling-free adaptive CORDIC rotator”  IEE Proceedings on Computers and Digital Techniques, 151, (6), 448-456.
Available at http://eprints.soton.ac.uk/263504/

Original foundation work on the CORDIC Algorithm:
Jack E. Volder, "The CORDIC Trigonometric Computing Technique," IRE Transactions on Electronic Computers, Vol. EC-8, pp.330-334, 1959.
John S. Walther, "A Unified Algorithm for Elementary Functions," Spring Joint Computer Conference Proceedings, Vol. 38, pp.379-385, 1971.
D. H. Daggett, "Decimal-Binary Conversion in CORDIC," IRE Transactions on electronic Computers, Vol. EC-8, No. 3, pp. 335-339, 1959.
M. A. Liccardo, "An Interconnect Processor with Emphasis on CORDIC Mode Operation," Masters Thesis EE Dept, University of California at Berkeley, September 1968.
J. E. Volder, "Binary Computation Algorithms for Coordinate Rotation and Function Generation," Convair Report, IAR-1 148, Aeroelectronics Group, June 1956.

en.wikipedia.org/wiki/CORDIC:  Matlab/Octave code for v=cordic(beta,n)


How to do sliding window multiplication in frequency domain

Is there a way to get my vectors Y1, Y2, Y3, etc... by computations in the frequency domain instead?

Yes. Rough outline: you take a full length FFT if the zero padded window. For each frame you can implement the time shift a multiplication with $e^jn\omega$ and then do a circular convolution and downsample. If you choose the window size well, you can do a down-sampled convolution instead.

And how do I know if it's faster/slower than just taking the inverse FFT of the frequency response and doing the steps described above?

It will almost always be slower. Convolution has computational complexity of $N^2$ whereas multiplication has only $N$. Saving a few short FFTs will not compensate of this.


Filtering difference of correlated measurements

I have a couple thoughts on your problem:

Regarding the handling of atmospherically-induced biases to the GPS receivers' reported positions, I don't think you need to do anything at all. If there is a truly common bias inherent in both positions, then it will simply cancel when you take the difference of the two positions. Thus, I would just treat each platform independently and subtract the time-aligned positions to yield your displacement vector estimate.

This is essentially the technique of differential GPS (DGPS), which was used successfully to combat the effects of the US government's policy of selective availability (SA) on the publically-available GPS carrier. SA was implemented by pseudorandomly dithering the timebase of the spacecraft; this would result in small position errors reported by GPS receivers.

This type of implementation meant that the SA error was highly correlated with the receiver position, such that two receivers in close proximity would experience approximately the same SA error. In the typical DGPS case, if one receiver's position is known precisely, one can solve for the SA error and resolve the second receiver's absolute position accurately, cancelling out the effects of SA. In 2000, the US government ended the use of selective availability, so this isn't much of an issue any more.

Your case is similar but slightly different, in that you don't know the absolute position of either receiver, but that doesn't seem to be what you care most about. Instead, you only care about the differential displacement between the two. Therefore, you can apply a similar technique; if there is a constant position bias inherent in both, then calculating the difference of the two positions will cancel it out, just as DGPS systems did for SA error.
It's not likely that you can improve the performance of the GPS receivers' reported positions by filtering them further. Dynamics tracking systems are a key part of any GPS receiver, and a lot of work has gone into their development. Especially if the receivers that you are using are of a recent vintage, I would hazard that any additional filtering that you would apply would be superfluous.

That's not to say that you would never want to process the reported positions further. There are some exceptions to my above statement:

If you have data from other instruments to fuse with the GPS readings, that can yield a gain. Sensors used for dead reckoning (e.g. gyroscopes or accelerometers) can yield measurements that can help to improve tracking of platform dynamics in periods where the GPS receiver outputs are unreliable or inaccurate. However, in time periods where the receiver measurements are assumed to be accurate, an optimal tracker would rely proportionally more upon the GPS measurements in lieu of the higher-variance dead reckoning sensors.
If you have some specific knowledge of your platform's dynamics that you know the GPS receiver isn't taking into account, then you may be able to obtain some gains. If you can specifically configure your receiver for maritime operation, you may be able to obtain better performance, as its tracker can make certain assumptions about your platform's position (e.g. zero altitude relative to sea level, no vertical velocity or acceleration). Failing that, it's possible that you could take the receivers' position outputs and apply some filtering to improve things, but your results will be very specific to the equipment that you use and whatever processing is applied by their internal trackers.


For a sine wave with a period $t_p$ the clock jitter or clock edge time error $t_j$ gives a phase error of $2 \pi t_j/t_p$ radians, where $2\pi$ represents a full period of a sine wave in terms of phase. The period $t_p$ of the sinusoid is an inverse of its frequency $f$, so the phase error can also be written as $2 \pi t_j f$. No matter if error is the peak error or the root mean square error, or whatever, as long as it is the same kind for both clock edge time error (jitter) and phase.


down sampling and a factor of 2 without loss of information

If you down-sample without filtering, you will fold spectrum above half the new sample rate together with spectrum below half the new sample rate.  If both those spectral bands (above and below) are sparse enough (big gaps in the spectral content) in the right places, then it is possible that when folded, there will be no overlap of the folded non-zero high band spectrum with the non-zero low band spectrum, thus no loss of information due to aliasing.  If you know where the spectral gaps were, you could unfold the data (with some effort).

In the more typical case where those two spectrum bands (upper and lower) are not sparse enough and thus do overlap, the two spectrums be mixed (like mixing paint) when downsampling, thus incur an informationally lossy process.


How to find the phase spectrum of a rectangular pulse? (Fourier Transform)

For the rectangular pulse without shifting, the phase spectrum its quite simple to obtain, since its Fourier Transform $F(\omega)$ is a real function. Real functions may only lay on the real axis, thus having $0$ or $\pm \pi$ phases, depending on the positive or negative sign of the function, respectively. Try visualizing the complex plane to realize this.

Now, if you're new to DSP, for $g(t)$ it might be tricky. First, take into account the scale change when dividing by 2 in $g(t)$. This expands the spectrum by a proportion of 2, i.e. $f(t/2)$  transforms as: $2~F(2\omega)$. That's why the $sinc$ function in $g(t)$ is twice as wide as in $f(t)$. Then the amplitude is multiplied by 3, that's why $G(\omega) = 6$ for $\omega = 0$.

Due to time shifting, you'd need to multiply $F(\omega)$ by a complex exponential. Recall the time-shifting Fourier transform theorem: $f(t - t_0)$ transforms as $F(\omega)~e^{-j\omega t_0}$.

The amplitude spectrum isn't modified by a time shift (since $|e^{-j\omega t_0}| = 1$), but the phase spectrum is added to $-\omega t_0$, which is the phase of the complex exponential (i.e., $\angle e^{-j\omega t_0} = -\omega t_0$).

Finally, you'd only need to add to the phase spectrum of $f(t)$ the linear function of $\omega$: $\angle H(\omega) = -4\omega$ (since $t_0 = 4$ in this case), which has a negative slope.

That's why you have in Fig. (d) the sum of a negative-slope linear function and a rectangular pulse train of $\pm \pi$ values.

This is basic DSP, I'd recommend to take a look to these theorems at Wikipedia for a very first glance: https://en.wikipedia.org/wiki/Fourier_transform#Properties_of_the_Fourier_transform


Why complex signal has no imaginary spectrum

The imaginary part of the spectrum corresponds to the odd part of the time-domain sequence. Since the given signal is even, i.e.,

$$x[n]=e^{jn\omega_0}=x^*[-n]=e^{-j(-n)\omega_0}=e^{jn\omega_0}$$

the imaginary part of the spectrum is zero, i.e., the spectrum is purely real-valued.

Note that for complex-valued signals, even in this context means $x[n]=x^*[-n]$, i.e., its real-part is even and its imaginary part is odd.

In sum, if a sequence $x[n]$ satisfies

$$x[n]=x^*[-n]$$

(i.e., the sequence is even and, consequently, its odd part is zero), then its DTFT is purely real-valued.


Why does the BER rate change randomly versus the SNR?

There are a few problems with your code:

You are using the index ii not the actual SNR in DB in this line: y= awgn(sFilt,ii,0);
You are not taking into account the filter delay in ySamp = yFilt(os:os:N*os);.  Your filter is filter = rcosine(1/T,os,'sqrt',rolloff); which will design a root-raised cosine FIR filter.  Your filter delay will be (length(filter)-1)/2.


I created a simple example with Scilab/XCOS. It is a simulation of lowpass filter with noise. It is observed by a kalman filter. You can also insert some uncertainties in the system model. With XCOS you can simulate the system. The results will be plotted automatically.
Here is the Downloadlink.

Little help with scilab:
Go to the directory with standard unix command

cd /home/workingdir/Kalman_Example

There you should execute getd() to load all functions (.sci-files) in the directory.
Then you can execute the modeling.sce  file, which is just an initializing script:

exec('modelling.sce')

Then typexcos simulation_kalm1.zcos into the console. A graphic window will open. Press the play button and the simulation will start.


Incoherence: Compressed Sensing (CS) vs Matrix Completion (MC)

Partial answer:

2) Matrix completion does "not rely on any notions of restricted isometry." Candes and Recht, 2008.


The output signal will still be normally distributed, but its power spectrum, i.e. its frequency content, will obviously be different from the input signal. If $S_X(\omega)$ is the power spectrum of the input signal, which is approximately flat, then the power spectrum of the output signal is

$$S_Y(\omega)=|H(\omega)|^2S_X(\omega)$$

where $H(\omega)$ is the frequency response of the filter.


Calculating information loss in a signal after smoothing?

Information loss and KL Divergence are functions that describe the behavior of random variables, in other words a set of random signals from some distribution. When you only have one signal, like you do in this case, calculating the information loss doesn't make sense. You can calculate mean squared error like the other comment suggests.


Generally speaking, a change of color space is essentially a nonlinear transform that stretches and distorts the coordinates. The transform is also continuous, except maybe on a singularity line or half-plane. It does not preserve the distances, but preserves the neighborhoods.

In the case of RGB > Lab, you can see it as a linear transform (change of basis, from RGB to XYZ) followed by a nonlinear transform of every coordinate independently (rescaling of XYZ to Lab). The topology is essentially preserved.

On the other hand, image segmentation algorithms do their best to delimit clusters of similar colors and separate the clusters from other clusters. So if you just deform the clusters, a good segmentation algorithm will adjust the boundaries and somehow undo the transform. If you repeat the experiment with XYZ or YUV coordinates, you should see no big difference either.

A stronger qualitative difference can be observed if the topology of the space is modified, for example by moving to the HSV space (Cartesian to cylindrical coordinates), of by performing some projection to lower the dimensionality such as XYZ > xyz. But beware of the physical interpretation of what you are doing.


Carrier Recovery, Frequency and Phase offset in (CP)FSK

If the two nodes run the same clock i would know perfectly where the signal starts from the end of the first chirp, after the fixed period of pause. In this case the signal i try to demodulate would have no phase shifting, correct?!

Wrong. The signal travels a distance. The distance it travels inherently means you get a phase shift – a full period ($2\pi$) for every wavelength (that's kind of the definition of wavelength, isn't it).

The signal gets reflected. Reflection often also leads to a phase shift.

The system with which you produce and measure the physical signal has a phase-shifting behaviour, too.

The only case where you'd not see a time shift is when all these effects add up to a multiple of $2\pi$ (because periodicity).

  If my clocks run differently, i will reconstruct a different signal, right? If yes, why is shifted in phase?

It's shifted in frequency. Frequency is the derivative of phase. That means that the phase shift is now a linear function of time.

  Can you point me out how to adjust, if possible, this offset or to know where to sample the signal?

To correct a statically phase-shifted signal, you'll need phase correction. To correct a frequency-shifted shifted signal, you'll need frequency correction.

In analog thinking, you typically employ phase locked loops (PLLs) of first order (static phase correction only) or higher order (frequency shift, or even frequency drift).

  Is this effect something like a Doppler?

The Doppler effect leads to a shift in frequency, too. But it's a different thing.


What is formula for quantization step size?

The actual quantization step size depends on the implementation of the quantizer, and would be by definition the input voltage change that would cause the digital output to increase by 1 bit.
You can find the specifics for the actual device you are using on the manufacturer datasheet, such as this example from ADI's AD9459 ADC, where the input range is from $\pm V_{ref}$ where $V_{ref}$ is an externally supplied reference voltage, and OR in the table is an out of range indicator:

The actual threshold down to 0.5 lsb levels is of no practical use since the self-noise within most ADCs typically exceed 1 lsb (otherwise the core could be used as a higher bit converter). Datasheets will also specific ENOB (Equivalent Number of Bits) which is dependent on sampling rate and input frequency. This captures all noise sources and provides a metric of the same number of bits the device provides compared to a "perfect" ADC. Note that this precision is maximized with an input signal that is typically 1 to 2 dB below actual full scale where clipping would occur (as we approach the actual full scale, spurious signals increase due to non-linearities which degrades overall dynamic range and therefore precision).
That said it is sufficient in most cases to simply use the total input range (up to the backed off full scale) and divide by $2^{b}$ where $b$ is the stated precision in ENOB.


Having difficulty checking for time invariance of discrete system

You did everything correctly. I think that part of your confusion comes from the way you chose the notation. Let me use $v_k[n]$ to denote the system's response to a delayed input $x[n-k]$. $v_k[n]$ is described by the following difference equation:

$$v_k[n]=3v_k^2[n-1]-nx[n-k]+4x[n-k-1]-2x[n-k+1]\tag{1}$$

On the other hand, the delayed response $y[n-k]$ to the original input $x[n]$ is described by

$$y[n-k]=3y^2[n-k-1]-(n-k)x[n-k]+4x[n-k-1]-2x[n-k+1]\tag{2}$$

The system would be time-invariant if $y[n-k]$ and $v_k[n]$ satisfied the same difference equation. From $(1)$ and $(2)$ we can see that that is not the case, hence the system is time-varying. You can see that the multiplicative term $n$ is responsible for making the system time-varying. In general, if there are time-dependent coefficients in a difference equation then the corresponding system is time-varying.


Just apply:

[B,A]=butter(3,2*f0/fs,'high');
y=filter(B,A,x);

where fs is your sampling frequency, f0=5 (your red band-width, as shown approximatedly in your spectrogram), x is your nx1 column vector data and y is your new filtered data.

A | B is the denominator | numerator of your highpass filter.

Hence the original picture:

spectrogram(y);

will show the 0-5Hz band in a pure blue band, a sign of deleting everything in there...

Best regards....


What algorithm can I use to "match" or "align" (translate and scale) two vectors of different size in C#?

The most similar method to your thought is finding a "linear regression". Hence, you can use a c# library or code to compute "linear regression" for the points of in A|B elements pairs.


Concept of Orthogonality and Orthonormality

Orthogonal means that the inner product is zero. For example, in the case of using dot product as your inner product, two perpendicular vectors are orthogonal. Orthonormal means these vectors have been normalized so that their length is 1.

Orthogonal vectors are useful for creating a basis for a space. This is because every point in the space can be represented as a (linear) combination of the vectors. So for example in 3D space, x=[0,0,1]  y=[0,1,0] and z=[1,0,0] form an orthornormal basis. No component of x can be represented with a component of the other vectors. This is because they are linearly independent

This is a different type of independence to statistical dependence however.

Correlation is a different notion. There are different ways of representing the correlation of two vectors (random variables X and Y).

$corr(X,Y) = \frac{cov(X,Y)}{\sigma_x \sigma_y} = \frac{E[X - \mu_x]E[Y - \mu_y]}{\sigma_x \sigma_y}$

Covariance is a measure of how two variables change together.

Another measure of correlation is cross-correlation. This is measure of the similarity of two functions (or vectors).

These concepts ares used in many different ways in image processing.

For example cross-correlation is used in template matching. When looking for a small image inside a much larger image, the small template image is 'slided' over the large image and the cross-correlation is computed for each position. Locations with a high cross-correlation are likely to contain the image we are looking for.

The concept of Linearly independent component vectors is used in principal component analysis (PCA). PCA takes a cloud of points in space and calculates a set of orthogonal basis vectors to represent them.

Independent component analysis (ICA) is used for separating mixed signals. ie: separating a speaker at a noisy cocktail party. ICA uses properties of the signal correlation to separate the signals.


The frequency resolution is independent of the sampling rate. Consider that in the extreme of an infinite sampling rate (continuous time), the Fourier Transform of a single FSK symbol at either 1070 or 1270 Hz carrier will be a Sinc function with the first null at $1/T$ away from the carrier where $T$ is the symbol duration. For a 300 bps data rate, the first null is offset from the carrier by 300 Hz. Clearly this is not orthogonal frequency shift keying since the other symbol is only 200 Hz away, but that doesn't mean the higher transmission data rate can be used, it just means the normalized peak to peak amplitude of the signal out of the frequency discriminator will be reduced due to the cross correlation between symbols.

The OP's sampling rate of 8 KHz is more than sufficient according to Nyquist. An optimum strategy would be to center two bandpass filters (Goertzel is a reasonable approach) at each of the tone frequencies and ensure the bandwidth of each filter is wide enough to allow for the bandwidth of the modulation rate to pass through, which is approximately the data rate, with care that the impulse response of any filtering doesn't cause inter-symbol interference.  The larger of the two bins in any given symbol period would be the most likely symbol that was transmitted.

Example 300 bps BFSK Demodulator

Below shows a complete non-coherent FSK demodulator implementation that works up to the higher 300 bps data rate and is very robust against carrier frequency offset (CFO); able to operate with offsets as high as the data rate. This implementation offers a nice balance of efficiency with performance.

The input modulated signal is frequency translated to baseband by multiplying with $e^{-j 2\pi f_c t}$ where $f_c$ is $1170$ Hz, midway between the FSK symbol frequencies. The creates a complex baseband I and Q output which is then decimated by 4, and the decimation filters efficiently remove the high frequency image of the modulated signal, and then provide a more efficient lower sampling rate of 2 MSps for the rest of the receiver processing. The decimation filters as half-band filters are linear phase filters with every other tap coefficient equal to zero, so very few actual multipliers are required. The code demonstrated here with 17 tap filters requires only 5 multipliers for each filter and provides the following performance in the plots given below, maximizing rejection of the larger high frequency noise component from the frequency translation stage.

The four-quadrant inverse tangent (ATAN2(I/Q) function is used to derive the instantaneous phase of the baseband signal, and then this phase is differenced (approximating the ideal differentiation which converts phase versus time to frequency).

The differenced phase is then unwrapped (so that phase doesn't jump suddenly between $+\pi$ to $-\pi$ when the $\pm \pi$ boundary is passed)  by detecting $|\Delta \phi|>\pi$, and adding or subtracting $\pi$ accordingly. The unwrap() function in MATLAB does this under the hood, but convenient here to implement directly since we are taking the difference between successive samples anyway.

The mean value of the unwrapped differenced phase will be proportional to any carrier offset frequency assuming equiprobable data with out an excessively long run of all 1s or all 0s. These factors, as well as the possible rate of change of the carrier offset in dynamic conditions and the acceptable overall delay in the receiver is then considered when determining the number of samples to average over, and then this average is subtracted from a delay-matched copy of the processed signal (matching the delay is only important if the demodulating the initial bits while the offset is being determined is desired, at the expense of overall latency).

Finally large overshoots that may occur at the frequency transitions (especially if the modulation does not transition smoothly from one frequency to the next) are removed through a tanh function, which is reasonably approximated by:

$$tanh(x) \approx \begin{cases}x,& |x|<1\\1,& x>1\\-1,&x<-1 \end{cases}$$

For larger SNR's this can simply be $sign(x)$, but using the mapping above will result in better performance in lower SNR conditions.

After any carrier offsets and larger transition spikes are removed the signal can be further filtered to optimize SNR performance prior to decision.  This filtering can be incorporated into the timing recovery, and a very efficient but high performance approach would be to resample with polyphase filters which can provide filtering and fine delay adjustment and operate together with a timing detector as part of a timing recovery loop (for example, see  Gardner Timing Recovery for Repeated Symbols and What is the difference of each frequency response of partial filters in a polyphase method).

Here is the result with a random data sequence followed by the demonstration MATLAB/Octave code

Note: it is very important to test with a random data pattern containing both rapid 10101 transitions as well as longer runs. Testing with a repeating 1010 pattern will confirm highest data rate but will also be the most behaved result. Below shows the equivalent result using the OP's modulated data with this demodulator to a repeating 1010 pattern:

Code:

Example modulated data  with worst case hard transitions between f1 and f2

fs = 8000;
f1 = 1070;
f2 = 1270;
rate = 300;
#create example modulated data
syms = [1 -1 1 1  -1  1 1 1 -1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1];
N = floor(length(syms)*fs/rate);
n =0:N-1;
t= n/fs;
data = ones(1, length(n));
sym =  1
for i = [1:length(n)]
    data(i) = syms(sym);
    if (i > floor(sym*fs/rate))
      disp(syms(sym))
      sym = sym+1;
    end
end
fc = (f1+f2)/2;
fd = fc-f1;
fco = 0;     # set to sim carrier freq offset
modx = cos(2*pi*(fc+fco-data.*fd).*t );

FSK Demodulator:

# modulated data at 8 KSps in array modx

t = [1:length(modx)]/fs;  # time vector for array

# Translate to baseband
demod = modx.*exp(-j*2*pi*(f1+f2)/2*t);

# Decimation filter
ord = 16;        #set to multiple of 4
k = (-ord/2:ord/2);
coeff = 0.5*sinc(k/2);
win = kaiser(ord+1,5);
coeff = coeff.*win';

# Half band decimations
demod2 = filter(coeff, 1, demod);
demod2= demod2(1:2:end);
demod4 = filter(coeff, 1, demod2);
demod4= demod4(1:2:end);

# ATAN2 and derivative filter
Q= imag(demod4);
I= real(demod4);
delta = filter([-1 1],1, (atan2(Q,I)));

# Unwrap phase, same as unrwap(atan2(Q,I))
delta2 =  mod(delta+pi,2*pi)-pi;

#subtract mean (CFO removal)
out = (delta2- mean(delta2));

# normalize and eliminate peaking at transitions
outlim= tanh(5*out/std(out));

Another approach that is even simpler is to use a delay and multiply frequency discriminator as as I had further detailed in this post: https://electronics.stackexchange.com/questions/293706/fsk-demodulation-using-dsp/293723#293723. I was not able to create an implementation anywhere close to the performance of the version above at the higher 300 bps rate since after the delay of 19 samples there are only 6 samples remaining within the symbol for the product. However an implementation is further detailed below with good results for the lower 110 bps rate.

The digitally implemented delay and multiply discriminator approach from the linked post is copied below:

Specific to this application for a real signal, the output of the low pass filter is proportional to the cosine of the phase difference at the input to the multiplier Therefore, the delay $T$ is set such that the average frequency 1170 Hz is an integer multiple of 90°, and for a maximum response $T = 1/(2 f_\Delta) = 2.5$ ms, which is 20 samples at 8KHz. At 1170 Hz this is 117°; reducing the delay by one sample decreases T by 125 us and moves the 1170 Hz crossing to $ -79.7°$ best centering the FSK signal in the demodulator within an integer sample delay. This is shown in the graphic below where the blue line represents the discriminator response for a 19 sample delay versus input frequency: When the input is at 1070 Hz the normalized discriminator output is close to -1, and when the input is at 1270 Hz the normalized discriminator output is close to +1.

Frequency Discriminator Response of Delay 19 Samples and Multiply

The multiplier produces the sum and the difference of the input frequency for a real signal. A moving average filter can be used as an efficient low pass filter rejecting the sum products and passing the difference (which is the demodulated signal through).

Resulting FSK Demodulator

Thus for comparison to the Goertzel, we have a complete non-coherent FSK demondulator for a real modulated FSK 1070 Hz/1270 Hz signal sampled at 8KHz, implemented with only a 19 tap delay, 1 real multipliers and 23 adders. The noise bandwidth of this filter is 350 Hz, so with additional filtering SNR can be improved further. Cascading this with a polyphase resampling filter to twice the symbol rate would provide additional filtering with minimum processing, and can be used as the timing adjustment in a timing recovery loop (such as the Gardner loop which operates at 2x samples per symbol).


Why can Quadrature Demodulation demod a Frequency Modulated Signal?

I don't think your formula is quite right.

The two complex samples are the locations (at two successive sampling instants)
of the tip of the rotating phasor that represents the analog signal.  The information that you need is the angle between phasor positions at these
two successive time instants.  If $S_{n+1} = r_{n+1}e^{j\theta_{n+1}}$ and
$S_n = r_ne^{j\theta_n}$, then you want the value of $\theta_{n+1}-\theta_n$.
Now, $S_{n+1}S_n^* = r_{n+1}r_ne^{j(\theta_{n+1}-\theta_n)}$ and if you
express these complex numbers in rectangular coordinates, then you can
get the desired angle as
$$\theta_{n+1}-\theta_n
= \arctan\left(\frac{\text{Im}(S_{n+1}S_n^*)}{\text{Re}(S_{n+1}S_n^*)}\right)
= \arctan\left(\frac{\text{Re}(S_{n})\text{Im}(S_{n+1})
- \text{Im}(S_n)\text{Re}(S_{n+1})}{\text{Re}(S_{n})\text{Re}(S_{n+1})
+ \text{Im}(S_{n})\text{Im}(S_{n+1})}\right)$$
though usually the four-quadrant arctangent function atan2(y,x) is used
instead of atan which returns an answer between $-\pi/2$ and $\pi/2$.

In digital communications applications such as DQPSK demodulation,
the real interest is not in the actual value of the angle, but in which
of the four intervals $(-\pi/4,\pi/4)$, $(\pi/4, 3\pi/4)$, $(3\pi/4, 5\pi/4)$,
and $(5\pi/4, 7\pi/4)$ the angle belongs -- in other words,
did the phase not change at all, or change by $\pi/2$, or by $\pi$, or
by $3\pi/2$ -- and this can be readily determined
by comparing the values of the numerator and denominator in
the argument of the arctangent above, thus saving computational
effort and speeding up the demodulation process. But the
more general formula gives the actual angle, and thus can be
used to demodulate a general FM signal as well.  The output
will be a sequence of samples of the original continuous-time
signal that modulated the carrier and thus formed the signal
that was transmitted.  D/A conversion will give the reconstructed
signal.


Real world application of signal sparsity?

Sparsity concept is extensively being used in computer vision and image processing. The Idea is that natural image can be pretty sparse when it is transformed to different bases. this bases can be predefined, e.g. FFT,DCT or can be learned from the image, e.g. sparse coding.

Here are a few well know example of algorithms that uses the sparsity assumption on images:

Jpeg image compression: Uses DCT transformation. assumes that high frequencies are low in natural image, so they are being removed in order to preserve space for compression. see more information on Wikepedia
Sparse coding: is a technique being used that assumes that an image has a lot of redundant information and because of that can be reconstructed from an overcomplete dictionary. each patch in the image should be sparse combination in the dictionary. see instroduction to sparse representation presentation
inverse problems in computer vision like recovering image embedded in noise and recovering sharpe image from blurred image are also using sparse assumption, e.g. wavelete shrinkage is being used for image denoising. See application to wavelet trasform in here


Hilbert Transform in C provides possibly strange results

The function in your question doesn't compute the Hilbert transform.

The Hilbert transform can be implemented either in the time domain by filtering the input signal with a filter kernel that approximates the ideal kernel of the Hilbert transform (there are approximations for FIR and IIR Hilbert transformers), or it can be implemented in the frequency domain.

An example of a frequency domain implementation is the Matlab/Octave function hilbert. But note that this function actually computes the analytic signal, so the Hilbert transform is the imaginary part of the result.

A simple direct frequency domain implementation of the Hilbert transform is shown in the Matlab/Octave script below. Note that it only uses real (I)FFTs.

% compute Hilbert transform of vector x

N = length(x);
X = fft(x);

if ~rem(N,2),       % N is even

    n = N/2;

    % set values at DC and Nyquist to zero
    X(1) = 0; X(n+1) = 0;

    % multiply with -1i * sign(w)
    X(2:n) = -1i * X(2:n);
    X(n+2:N) = 1i * X(n+2:N);

else                % N is odd

    n = (N+1)/2;

    % set value at DC to zero
    X(1) = 0;

    % multiply with -1i * sign(w)
    X(2:n) = -1i * X(2:n);
    X(n+1:N) = 1i * X(n+1:N);

endif

y = real(ifft(X));


Antennas don't just work for a single frequency and for nothing else. They typically work for a larger range of frequencies. We call the total amount of frequencies their bandwidth.

So although FSK uses a set of carrier frequencies, all of these are typically well-received by a single antenna.

Yes, bandwidth of antennas is limited, but so is the bandwidth of the signals you want to send over the air or can even generate and process.

"Radio Stations": you probably mean FM audio broadcasts? Those are FM, not FSK; the same applies, however: though FM uses a whole bandwidth to transport a single programme, an antenna can easily pick up this whole bandwidth. In fact, normal "FM Radio" antennas capture the whole ca. 87 MHz – 108 MHz bandwidth reasonably well at all.

You seem pretty confused by a set of basic terms. I don't really know your background, but getting a clear introduction to what a signal, frequencies, bandwidth, and radio are might be a good idea before getting into concepts like FSK. In fact, I'm pretty confused by the fact that you already asked about LDPCs and now about these basic terms of radio transmission – I'd typically expect people to learn about what an antenna does and how that works with signal bandwidths a solid year before they can get to learn about LDPCs – which demands good understanding of what a radio channel is. So my intuition is that you're from a math or computer science background, and it came easier to you to understand algebraic coders than  RF channels (which would explain why you didn't immediately understand, though you showed great insight there, that deletion channels aren't that often used to model wireless communication).

Now, aside from getting a nice textbook (or the ARRL handbook, for a start) that explains a bit of RF technology basis, it might be a good idea for future questions to explain a bit of the background you have. I'm really puzzled by the "order" of your questions!


The $\sigma$ decides the scale of objects being simplified. This is explained here:

  The size of the Gaussian filter: the smoothing filter used in the first stage directly affects the results of the Canny algorithm. Smaller filters cause less blurring, and allow detection of small, sharp lines. A larger filter causes more blurring, smearing out the value of a given pixel over a larger area of the image. Larger blurring radii are more useful for detecting larger, smoother edges – for instance, the edge of a rainbow.

Thus the support size of the Gaussian function decides the level of details that are simplified. Simplest way is to test with binary squares/circles.


DFT of real sinusoids - why sum over -$N/2$ to $N/2-1$ as opposed to $0$ to $N-1$?

The DFT (Discrete Fourier Transform) of an $N$-point sequence $\left\{x[n]\right\}$ is given as:
$$X[k]=\sum_{n=0}^{N-1} x[n] e^{-j2\pi kn/N}, \quad k=0, 1, \ldots, N-1\tag{1}$$

If we assume the sequence $\left\{x[n]\right\}$ and its DFT are complex, they can be expressed as follows (with the subscript $r$ and $i$ for real and imaginary part respectively):

\begin{align}
x[n] &= x_r[n]+jx_i[n], \quad 0\leq n\leq N-1\\
X[k] &= X_r[k]+jX_i[k], \quad 0\leq k\leq N-1
\end{align}

With a bit of algebra and using equation $(1)$, you get:

\begin{align}
X_r[k] &= \sum_{n=0}^{N-1}\left[ x_r[n]\cos\left(\frac{2\pi kn}{N}
\right)+x_i[n]\sin\left(\frac{2\pi kn}{N}
\right)\right]\\
X_i[k] &= -\sum_{n=0}^{N-1}\left[ x_r[n]\sin\left(\frac{2\pi kn}{N}
\right)-x_i[n]\cos\left(\frac{2\pi kn}{N}
\right)\right]
\end{align}
If the sequence $\left\{x[n]\right\}$ is real-valued, this means $ x[n]=x_r[n]$ as $x_i[n]=0$. The DFT will then be expressed as:

$$X[k]=\sum_{n=0}^{N-1} x[n]\left[\cos\left(\frac{2\pi kn}{N}
\right)-j\sin\left(\frac{2\pi kn}{N}
\right)\right]\tag{2}$$

And it can be shown from equation $(1)$ and $(2)$ that for these real-valued signals:

$$X(N-k)=X^*(k)=X(-k)\tag{3}$$

From the above results, it suffices to compute the DFT for half of the values of $k$ to get the full characteristic of the spectrum. From the $N$ values you compute less (roughly half) $k$ values. The boundaries of the summation are still the same, it still goes from $n=0$ to $N-1$ but your $k$'s shouldn't go up to $N-1$ as it will be a waste of resources.
$$X[k]=\sum_{n=0}^{N-1}y[n] e^{-j2\pi kn/N}$$

For $N$ even, $ 0\leq k \leq N/2$ suffices. With the results from equation $(3)$ you can get the remaining values from $N/2 + 1$ to $N-1$ or from $-N/2 + 1$ to $-1$.
You have $ 0\leq n \leq N$. So, $y[-N/2]$ does not exist and $y[N/2]$ is the value of your sequence $\left\{y[n]\right\}$ at $n=N/2$ (if $N$ is even).
However, for a real even sequence $\left\{y[n]\right\}$, the DFT $X[-N/2]$ equals $X[N/2]$ and this is the value at Nyquist.
Also note from equation $(3)$ that if $N$ is even, $X[0]$ and $X[N/2]$ must be real, and if $N$ is odd $X[0]$ is still real but you don't have a frequency bin at exactly $N/2$. Equation $(3)$ applies for real-valued sequences, what has been concluded from that equation doesn't apply for strictly complex-valued sequences. The advantage is that you've saved computations by computing only half of what you need to compute.


Is there a version of Welch's method that doesn't look for power?

Since an FFT is a linear operator, adding up the complex results of a sequence of FFTs of short windows is the same as doing a single short FFT on the the vector addition of all those short windows.

Note that for signals that are exactly integer periodic in the FFT width (sequential 0% overlapped windows), the vector addition will constructively interfere.  But for signals that are not exact integer periodic in the FFT width, there can be significant destructive interference.  So the result will be highly scalloped, with a narrow high gain bandpass around each FFT bin center, and potentially a lot of signal loss between bin centers.  Sometimes this is useful, sometimes not, but a very different result from a single longer FFT over the entire window set.


Instantaneous frequency vs fourier frequency

the frequency of a sinusoid is normally thought to be undefined when there is no sinusoid (i.e. when it is turned "OFF").

now, there are at least a couple of different methods of turning it OFF.  the most obvious is to effectively multiply the waveform output by zero to disconnect the sinusoid from the output.

however if the output of your signal generator is always connected (so there is no multiplication by zero) to the look-up table (LUT) that has the waveform values in it (and let's say that those waveform values follow a sinusoidal function), then another way to stop the output is to stop advancing position in the LUT.  that would be like setting the instantaneous frequency to zero.  but, unless you're watching out for this, it could get stuck on a DC value that is non-zero.  so there are a lot of functions, besides the zero function, that have a frequency of zero.  any DC function has instantaneous frequency of zero.


Is there an efficient algorithm to extract low frequency signal from audio with high frequency carrier?

If you know your carrier is always exactly 16kHz (within a few percent or so), and you have sampled the audio at 48kHz, you can try a simple 3-tap FIR filter approach:

$$Smoothed(k) = \frac{1}{3}( Sig(k) + Sig(k-1) + Sig(k-2) )$$

This assumes your original signal is called $Sig()$, and it computes a new signal, $Smoothed()$.

The method exploits the fact that three sample periods equates to one entire cycle of a  16kHz sine wave, so an 16kHz component in the input signal will be cancelled (since summing all the samples of a full sine wave cycle will result in zero).

The resulting $Smoothed()$ signal should look like you original plots, with the furry 16kHz components gone.

If your sample rate is not 48kHz (let's call it $F_s$), or if your interfering sine wave is at some other frequency (let's call it $F_i$), you can craft the 3-tap filter as follows:

$$Smoothed(k) = \frac{1}{2-2\cos(\omega)}( Sig(k) - 2\cos(\omega) Sig(k-1) + Sig(k-2) )$$

where $\omega = 2\pi\frac{F_i}{F_s}$ (note that, if $F_i=16000$ and $F_s=48000$, we get $\omega=\frac{2\pi}{3}$ and hence $\cos(\omega)=-\frac{1}{2}$. If you substitute this into the second equation, you will simply get the first equation again.

This method is simply creating a $2^{nd}$ order FIR filter with zeros at $\pm F_i$, and with the gain at DC normalised to unity.


How to measure correlation between 1 bit samples and theoretical original waveform?

I've sampled it at a 1 bit resolution (as 8 bit oscilloscope readings & 0b1)

No. Taking the LSB is NOT 1 bit quantization. A one bit quantizers are not trivial to design. For a zero mean signal that's reasonably symmetrically distributed over positive and negative values, you can just take the sign bit. For signals that are asymmetrically distributed around their mean, designing a decent 1 bit quantizer takes a bit of work.
In your example you should quantize so that you get a "1" of the signal is above the mean and and "0" if it's below the mean. If you do this, you will see that both quantization noise and quantized signal are highly correlated with the original signal. Both will be periodic with the same period.
If you want to study quantization noise, I recommend starting with higher bit levels, i.e. where the amplitude of the signal is much larger than the quantization steps. This makes the choice of quantizer very simple and in most case you can assume that the quantization noise is uniformly distributed. Otherwise, you need to tackle the subject of "optimum quantizer" which requires some non-trivial math.


Object-Recognition using SIFT: Switching the order of images creates two different results

ARGH. I feel ashamed. How about a cross validation? That's the solution. Now there are almost exact matches ... almost. =) Thank you, Libor. The second comment gave me a hint. =)


Which anti aliasing filter algorithm for efficient audio resampling?

For a first cut, I would consider a windowed Sinc interpolation kernel for both resamplers.  Very narrow window with a small polyphase table for the fast one.  Much wider window with more phases in a larger, possibly interpolated lookup table for the slower higher quality filtered interpolator/resampler.  Lots of parameters (including choice of window) that can be adjusted for both quality and performance to meet widely varying requirements.


Frequency response of a (rectangular) integrator

Those are good questions! Here's my "two cents." Hilmar is correct, the integrator is not stable because there's a pole on the z-plane's unit circle. But I'll bet a pint of beer that the region of convergence is outside the pole (and therefore outside the unit circle) because the integrator's impulse response is causal.

Answer to question# 1: Matlab's freqz() command computes the freq response at zero Hz by evaluating a ratio of polynomials at various frequencies. And that ratio has a zero-valued denominator at zero Hz. So freqz() gives a freq mag response of infinity at zero Hz.  The DFT of an impulse response that's a finite-length sequence of ones, evaluated at zero Hz, is merely the sum of the impulse response samples.  The two approaches result in freq magnitude responses that are zero-valued at all frequencies except at zero Hz.  Both approaches give similar results.

Answer to question# 2: You could, if you wish, call an integrator a "filter." But it's unstable, and unstable filters are of no value in trying to filter real-world signals.


What is the difference between undersampling and oversampling in analog to digital conversion ?

Good question as you can actually undersample and oversample at the same time! See my "DSP Puzzle" question on that specifically here:
How do you simultaneously undersample and oversample?

To best explain undersampling and oversampling, it is worthwhile understanding the concept of "Nyquist Zones" first.  This was explained in detail recently at this post:

Higher order harmonics during sampling

From which I am copying some graphics for further explanation here.

Sampling in the 1st Nyquist Zone, as in the graphic below, occurs when the real analog signal is located anywhere in the frequency range of DC to half of the sampling rate. Sampling at twice the highest frequency (and actually twice the highest bandwidth, as we will see more clearly when describing "under-sampling") is the minimum necessary condition to capture all information within that bandwidth according to the Nyquist-Shannon theorem. In practice to allow for achievable anti-aliasing filtering in the analog we need to sample a little higher than this, typically 20 to 30%.

The graphic below shows the related processing spectrums when sampling an analog input signal that is in the "first Nyquist Zone".  The colors are related to other discussion that is detailed in the link I gave above but not important here.

Oversampling

The above introduced standard "Sampling".  "Over-Sampling" is when we sample higher than the condition set by Nyquist, and typically much higher for us to be concerned with the possible benefit of doing that.

The main thing that happens when we over-sample is the quantization noise from the conversion gets spread over a much wider frequency range, thus in most cases reducing the noise from quantization in our frequency band of interest and therefore acting as a converter with a higher number of bits!  I have explained this further at this link so will refer to that for further detail on over-sampling:
What are advantages of having higher sampling rate of a signal?

Under-sampling

An analog signal can be higher in frequency than the first Nyquist Zone and still be sampled as long as the analog bandwidths of the system (including the A/D converter itself) are wide enough to pass the signal.  Such is the case in the graphic below copied for the link I first gave above that details this further. Here we are sampling a signal that is in the 3rd Nyquist Zone, and the resulting digital signal is indistinguishable from a similar signal that was in the 1st Nyquist Zone. So assuming a wide enough analog bandwidth, we can select any Nyquist zone to sample through the design of a anti-alias filter prior to the A/D converter. For standard and over-sampling conditions we would use a low-pass filter. And for under-sampling conditions we would use a band-pass filter. With no filter all the Nyquist zones would be sampled, thus we would have significant degradation due to noise and interference folding (aliasing). This is specifically the important reason to always have an anti-alias filter prior to the A/D converter!  "IF-Sampling" and "Direct-Sampling" are implementations of under-sampling. IF-Sampling is when we sample the IF frequency in a super-heterodyne receiver architecture, and Direct-Sampling is when we have an A/D converter with a high enough input analog BW (or a low enough input RF signal) that we can sample the RF signal directly.

Over-sampling implies having many more samples than the highest frequency of interest, and under-sampling implies we are down-converting the bandwidth of interest with a higher harmonic of the sampling clock (effectively). Thus as given in the first link, we can simultaneously under-sample and over-sample such as the following example using GPS:

Consumer GPS approximately occupies 2 MHz of BW at a carrier of 1575.42 GHz.  With an A/D converter that had an analog input bandwidth that exceeded this carrier frequency, we could directly sample the amplified and filtered signal from a GPS antenna with 100 MHz sampling rate (for example), which would result in a (spectrally inverted) copy of the GPS signal at 16* 100 MHz -1575.42 MHz = 24.58 MHz; well positioned for a digital IF in our 100 MHz sampling space. This would be digitally down-converted to baseband and then filtered where we would benefit from the advantages of over-sampling at 100 MHz the GPS signal that occupies significantly less bandwidth.  Thus we under-sampled by down-coverting with the sampling process, and over-sampled by choosing a sampling rate much higher than the minimum necessary.


Implementing a Matlab filter to cancel cross talk between two sensors

Try using other Blind source  separation(BSS) algorithms,even Kernel BSS could be applied to remove if it is nonlinearly mixed.


Object detection versus object recognition

You kind of answered your own question.

Object Recognition: which object is depicted in the image?

input: an image containing unknown object(s)

Possibly, the position of the object can be marked in the input, or the input might be only a clear image of (not-occluded) object.
output: position(s) and label(s) (names) of the objects in the image

The positions of objects are either acquired form the input, or determined based on the input image.

When labeling objects, there is usually a set of categories/labels which the system "knows" and between which the system can differentiate (e.g. object is either dog, car, horse, cow or bird).

Object detection: where is this object in the image?

input: a clear image of an object, or some kind of model of an object (e.g. duck) and an image (possibly) containing the object of interest
output: position, or a bounding box of the input object if it exists in the image (e.g. the duck is in the upper left corner of the image)


Inverse z transform - contour integration

In order to be able to do this, you need to know Cauchy's residue theorem, from which you get

$$\frac{1}{2\pi j}\oint_Cf(z)dz=\sum_kR_k\tag{1}$$

where $R_k$ are the residues of the analytic function $f(z)$ at its poles $p_k$, which must lie inside the closed contour $C$. The residue at pole $p_k$ is given by
$$R_k=\lim_{z\rightarrow p_k}(z-p_k)f(z)$$

In the case of the inverse $\mathcal{Z}$-transform, we have $f(z)=X(z)z^{n-1}$. Since in your example the region of convergence is given as $|z|>2/3$, we know that we need to find a causal sequence $x[n]$, so we only need to consider values $n\ge 0$. For $n>0$, there is only one single pole at $z=2/3$, and the corresponding residue is

$$R_0=\lim_{z\rightarrow 2/3}(z-2/3)X(z)z^{n-1}=-\frac{1}{3}\left(\frac{2}{3}\right)^{n-1}\tag{2}$$

For $n=0$ we get another pole at $z=0$ from the term $z^{n-1}$, and the corresponding residue is

$$R_1=\lim_{z\rightarrow 0}zX(z)z^{-1}=\lim_{z\rightarrow 0}X(z)=\frac12\tag{3}$$

According to (1), the value of $x[0]$ is given by the sum of the residues (2) and (3):

$$x[0]=-\frac13\left(\frac23\right)^{-1}+\frac12=0$$

For $n>0$ we only need to consider the residue $R_0$, so the inverse $\mathcal{Z}$-transform is given by

$$x[n]=-\frac13\left(\frac23\right)^{n-1}u[n-1]$$

Note that this is just meant as an exercise. In practice this example is a basic transformation, which you either know or look up in a table.