What is meant by *sampling* in terms of the *sampling theorem*?

I think, with respect to sampled function, that i disagree a bit with MBaz.  perhaps i misunderstand what he/she says regarding needing an infinite number of coefficients or samples for this special case of sampled periodic functions.  we had a little discussion about this at comp.dsp, and i'm gonna make use of $\LaTeX$ to spell out the same points.

of course, whether $x(t)$ is periodic or not, if $x(t)$ is real and is sufficiently bandlimited (there are no frequency components as high or higher than $\frac{f_\text{s}}{2} = \frac{1}{2T_\text{s}}$), then the samples:

$$ x[n] \triangleq x(n T_\text{s}) $$

are sufficient to completely represent the continuous-time $x(t)$.  if $x(t)$ never repeats, than an infinite number of discrete $x[n]$ are necessary to represent $x(t)$ for all $t$.

but if $x(t)$ is periodic,

$$ x[n+N]=x[n] \quad\quad \forall n,N \in \mathbb{Z} $$
and
$$ x(t+T)=x(t) \quad\quad \forall t $$

then samples existing over the span of one period are sufficient to represent $x[n]$ and $x(t)$.  in order for the $x[n]$ to be periodic in the sampled domain:

$$ \begin{align}
x[n+N] & = x\left((n+N)T_\text{s} \right) \\
& = x\left(nT_\text{s}+NT_\text{s} \right) \\
& = x\left(nT_\text{s}+T \right) \\
\end{align} $$
because $ x(t+T) = x(t) $.

which means the obvious, your function period $T$ has to be the same as $N$ times your sampling period $T_\text{s}$.

$$ T = N T_\text{s} $$

now, what we know about this sampled periodic function is that $N$ samples of $x[n]$ are sufficient to tell us all about $x[n]$, and since $x(t)$ is bandlimited, $x[n]$ and the $N$ samples that fully define it, are sufficient to fully describe $x(t)$.

if $x(t)$ is sufficiently bandlimited (as above), then

$$ \begin{align}
x(t) & = \sum\limits_{n=\infty}^{+\infty} x[n] \operatorname{sinc}\left(\frac{t-nT_\text{s}}{T_\text{s}}\right) \\
& = \sum\limits_{n=\infty}^{+\infty} x[n] \frac{\sin\left(\pi\frac{t-nT_\text{s}}{T_\text{s}}\right)}{\pi\frac{t-nT_\text{s}}{T_\text{s}}} \\
& = \sum\limits_{m=\infty}^{+\infty} \sum\limits_{n=0}^{N-1} x[n+mN] \frac{\sin\left(\pi\frac{t-(n+mN)T_\text{s}}{T_\text{s}}\right)}{\pi\frac{t-(n+mN)T_\text{s}}{T_\text{s}}} \\
& = \sum\limits_{n=0}^{N-1} \sum\limits_{m=\infty}^{+\infty} x[n] \frac{(-1)^{mN}\sin\left(\pi\frac{t-nT_\text{s}}{T_\text{s}}\right)}{\pi\frac{t-(n+mN)T_\text{s}}{T_\text{s}}} \\
& = \sum\limits_{n=0}^{N-1} x[n] \sin\left(\pi\frac{t-nT_\text{s}}{T_\text{s}}\right) \sum\limits_{m=\infty}^{+\infty} \frac{\frac{T_\text{s}}{\pi}(-1)^{mN}}{t-(n+mN)T_\text{s}} \\
\end{align} $$

$$ x(t) =
\begin{cases}
\sum\limits_{n=0}^{N-1} x[n] \frac{\sin\left(\pi\frac{t-nT_\text{s}}{T_\text{s}}\right)}{N\tan\left(\frac{\pi}{N}\frac{t-nT_\text{s}}{T_\text{s}}\right)}, & \text{if }N\text{ is even} \\
\sum\limits_{n=0}^{N-1} x[n] \frac{\sin\left(\pi\frac{t-nT_\text{s}}{T_\text{s}}\right)}{N\sin\left(\frac{\pi}{N}\frac{t-nT_\text{s}}{T_\text{s}}\right)}, & \text{if }N\text{ is odd}
\end{cases} $$

proving the latter takes a little bit.  you might recognize the $N$ odd case as the Dirichlet_kernel.  the $N$ even case looks a teeny bit different.  but both show exactly how the $N$ samples that fully define the sampled and periodic $x(t)$ are combined to get $x(t)$.

now, since $x(t)$ is also periodic with period $T$, then

$$ \begin{align}
x(t) & = x\left(t+T \right) \\
& = x\left(t+N T_\text{s} \right) \\
& = \sum\limits_{k=-\infty}^{+\infty} X[k] e^{i 2 \pi \frac{k}{T} t} \\
& = \sum\limits_{k=-\infty}^{+\infty} X[k] e^{i 2 \pi \frac{k}{N T_\text{s}} t} \\
& = \sum\limits_{k=-\lfloor \frac{N}{2} \rfloor}^{+\lfloor \frac{N}{2} \rfloor} X[k] e^{i 2 \pi \frac{k}{N} \frac{t}{T_\text{s}}}
\end{align} $$

where $\lfloor\cdot\rfloor$ is the $\operatorname{floor}(\cdot)$ operator and

$$ \begin{align}
X[k] & = \frac{1}{N} \sum\limits_{n=0}^{N-1} x[n] e^{-i 2 \pi \frac{nk}{N}} \\
& = \mathcal{DFT}\{x[n]\}
\end{align} $$

there's actually something to fudge (a factor of $\frac{1}{2}$) about $X\left[\frac{N}{2}\right]$ for the $N$ even case:

$$ \begin{align}
X\left[-\frac{N}{2}\right] = X\left[\frac{N}{2}\right] & = \frac{1}{2} \ \frac{1}{N} \sum\limits_{n=0}^{N-1} x[n] e^{-i 2 \pi n\frac{n(N/2)}{N}} \\
& = \frac{1}{2N} \sum\limits_{n=0}^{N-1} x[n] (-1)^n \\
\end{align} $$

note that:

$$ \begin{align}
x[n] = x(n T_\text{s}) & = \sum\limits_{k=-\lfloor \frac{N}{2} \rfloor}^{+\lfloor \frac{N}{2} \rfloor} X[k] e^{i 2 \pi \frac{k}{N} \frac{t}{T_\text{s}}}\bigg|_{t=n T_\text{s}} \\
& = \sum\limits_{k=-\lfloor \frac{N}{2} \rfloor}^{+\lfloor \frac{N}{2} \rfloor} X[k] e^{i 2 \pi \frac{k}{N}n} \\
& = \sum\limits_{k=0}^{N-1} X[k] e^{i 2 \pi \frac{nk}{N}} \\
& = \mathcal{iDFT}\{X[k]\}
\end{align} $$

if you deal with that $\frac{1}{2}$ fudging for $X\left[\frac{N}{2}\right]$ for the $N$ even case.  this is because $X[k+N]=X[k]$ for all $k$.

a periodic continuous-time function can be described with a countable set of Fourier coefficients.  a bandlimited periodic continuous function can be described with a finite set of $N$ Fourier coefficients, just as well as it can be described with a finite set of $N$ samples.


Consider an ideal low pass filter $H(\omega)$, and the input to this filter is the periodic square wave $x(t)$. Find the output $y(t)$

So you actually solve(d) this problem by computing the CTFS of the input $x(t)$:
$$ C_k = 5 \frac{ 1 - e^{ - j k \pi } }{ j k \pi } = \begin{cases} 0 &, \text{ for k even , } k \neq 0 \\ 10/jk\pi &, \text{ for k odd} \end{cases} $$ which gives the coefficient associated with the complex exponential of frequency $\omega_k = k \omega_0$ for the fundamental frequency $\omega_0 = \frac{ 2 \pi }{T_0} = \frac{ 2 \pi }{2}  = \pi$ , as the period of the input $x(t)$ is $2$ seconds.

Since the ideal lowpass filter passes those frequencies below its cutoff frequency $\omega_c = 4 \pi$, then we would only take those coefficiencs whose frequencies fall below $\omega_c$ yielding those values of $k = 0,\pm 1, \pm 3$ . Those indices higher than $\pm 5$ are filtered out by the ideal lowpass filter.

Therefore the output $y(t)$ is:
$$ y(t) = 5 + \frac{10}{j\pi} \left( e^{j\pi t} - e^{-j\pi t} \right) + \frac{10}{j3\pi} \left( e^{j3\pi t} - e^{-j3\pi t} \right) $$ which simplifies to:
$$ y(t) = 5 + \frac{20}{\pi} \sin(\pi t) + \frac{20}{3\pi} \sin(3\pi t) $$


This answer is a response to a comment by the OP on on yoda's answer.

Suppose that $h(t)$, the impulse response of a continuous-time
linear time-invariant system, has the property that
$$\int_{-\infty}^{\infty} |h(t)| \mathrm dt = M$$ for
some finite number $M$.  Then, for each and every
bounded input $x(t)$, the output $y(t)$ is bounded also.
If $|x(t)| \leq \hat{M}$ for all $t$ where $\hat{M}$
is some finite number, then $|y(t)| \leq \hat{M}M$ for all $t$
where $\hat{M}M$ is also a finite number.
The proof is straightforward.
$$\begin{align*}
|y(t)| &= \left |\int_{-\infty}^\infty h(\tau)x(t - \tau)\mathrm d\tau\right |\\
&\leq \int_{-\infty}^\infty |h(\tau)x(t - \tau)|\mathrm d\tau\\
&\leq \int_{-\infty}^\infty |h(\tau)|\cdot|x(t - \tau)|\mathrm d\tau\\
&\leq \hat{M}\int_{-\infty}^\infty |h(\tau)|\mathrm d\tau\\
&= \hat{M}M.
\end{align*}$$
In other words, $y(t)$ is bounded whenever $x(t)$ is bounded.

  Thus, the condition
  $\displaystyle\int_{-\infty}^{\infty} |h(t)| \mathrm dt < \infty$
  is sufficient for BIBO-stability.

  The condition  $\displaystyle\int_{-\infty}^{\infty} |h(t)| \mathrm dt < \infty$
  is also necessary for BIBO-stability.

Assume that every bounded input
produces a bounded output.  Now consider the input
$x(t) = \text{sgn}(h(-t)) ~\forall~ t$.  This is clearly bounded,
($|x(t)| \leq 1$ for all $t$), and at $t=0$, it produces output
$$\begin{align*}
y(0) &= \int_{-\infty}^\infty h(0-\tau)x(-\tau)\mathrm d\tau\\
&= \int_{-\infty}^\infty h(-\tau)\text{sgn}(h(-\tau))\mathrm d\tau
&= \int_{-\infty}^\infty |h(-\tau)|\mathrm d\tau\\
&= \int_{-\infty}^\infty |h(t)|\mathrm dt.
\end{align*}$$
Our assumption that the system is BIBO stable means that $y(0)$ is
necessarily finite, that is,
$$\int_{-\infty}^{\infty} |h(t)| \mathrm dt < \infty$$

The proof for discrete-time systems is similar with the obvious
change that all the integrals are replaced by sums.

Ideal LPFs are not BIBO-stable
systems because the impulse response is not absolutely integrable,
as stated in the answer by yoda.  But his answer does not really answer the question

  Can anyone give me a proof that ideal LPF can indeed be BIBO unstable?

A specific example of a bounded input signal that produces an unbounded output
from an ideal LPF (and thus proves that the system is not BIBO-stable)
can be constructed as outlined above (see also my comment on the main question).


If I build a low pass gaussian filter in the frequency domain should I also normalize it with the sum of the filter coefficients?

It is consistent to apply a normalization depending on what you really want to measure. Either on Fourier transforms, or low-pass linear filters, one customary sees three types of normalizations:

No normalization,
Normalization for a constant input,
Normalization in energy.

To be brief,  for a 3-point uniform coefficient filter, that is a low-pass filter, the impulse response $h$ can be:

No normalization: $h={1,1,1}$,
Normalization for constant output (mean-preserving): $h={\frac{1}{3},\frac{1}{3},\frac{1}{3}}$,
Normalization in energy: $h={\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}}}$.

Generally, the normalization depends on the invariance one wants to preserve. For low-pass filters, one often preserves mean (i.e. a constant image remains constant after filtering) or energy. For high-pass, some preserves energy, or estimated slope (or higher-order derivatives). Very often, the sum of coefficients of a high-pass filter is zero (a constant signal is zeroe out), thus it cannot be used for normalization. Yet, many other choices are possible.

So if you want to preserve the average value of the data, in time you normalize so that the coefficients sum to one. In frequency, you can remember that the amplitude of the $0$th frequency obtained from the Fourier transform is proportional to the sum of coefficients in time ("DC Term" in An Intuitive Explanation of Fourier Theory). So if it is non-zero, you can divise your whole Fourier spectrum by this value to get the same normalization.

The reason why people sometimes do not normalize the intensities with image filtering is that these values are often displayed and compared on a full range in the colormap (in Matlab, with function imagesc). Consequently, a scalar factor on the image is not considered harmful.


What is the difference between linear and non-linear filters?

To understand the linearity property more easily.Let us consider the above diagram,here we have 2 sequences namely Xn and Yn. when we add both the sequence we get Xn+Yn whose amplitude value are represented with blue colour. when any system which satisfy this condition then it is called linear.
In case of mean filter, mean value for sequence Xn is 1+1+3/3=5/3.mean value for sequence Yn is 1+2+0/3=1.mean value for Xn+Yn is 2+3+3/3=8/3.

$ mean(Xn)+mean(Yn)= mean(Xn+Yn), \\ \\
 5/3    +    1   = 8/3 $

hence we called mean filter as linear filter.
In case of median filter, if we calculate median value for sequence Xn, we get 1 (arrange the sequence in ascending order and then find the middle value). similarly median value of sequence Yn is 1 .median value of sequence Xn+Yn is 3.

$median(Xn)+median(Yn)  \neq  median(Xn+Yn).\\ \\
1 + 1  \neq    3 $

hence we call median filter as non-linear filter


Summary
You can have only one (or no) primary sensor which is used in the prediction process. All other sensors are secondary and only used for corrections. Each correction takes  the current state and a measurement from one of the secondary sensors and outputs the corrected state. You can run corrections as often as you want, typically each time one of the secondary sensors produces a new measurement. There is no need to run corrections at a fixed interval.
The new state is predicted from the model matrices F and B and from the old state x and primary sensor data (or control input) u.
$$x_{k+1} = Fx_k + Bu_k$$
The corrected state is calculated from the kalman gain matrix K, observation matrix H and secondary sensor data z. Each secondary sensor has its own observation matrix.
$$x_{k,corrected} = x_k + K_k(z_k-H_k x_k) $$
Detailed version
You have to designate one primary sensor each for orientation and position. Alternatively, you can go without primary sensors if system dynamics allow. The primary sensors should be extremely reliable, accurate in the short term, and support a high update rate. Long term drift is not a concern here.
In a car navigation system, the primary sensors are usually steering angle and odometer, the secondary sensor is GPS. In an attitude and heading indicator for aircraft, the primary sensor is a three-axis gyro and secondary sensors are accelerometer and magnetometer.
There can be an arbitrary number of secondary sensors. Their purpose is to provide long term accuracy, which the primary sensors lack. The secondary sensors have other shortcomings, like unreliability (no GPS in tunnels), poor short term accuracy (apparent gravity vector disturbed by g-forces from maneuvering), or incomplete information (gravity vector provides pitch and bank, but not heading). If you have a sensor that can take both roles - a sensor that is reliable, provides enough information to uniquely deduce orientation or position, and is accurate both in the short and long term, then you don't need a Kalman filter in the first place.
The prediction stage of the Kalman filter runs at a fixed and high rate and continuously updates the current estimates for position and orientation based on the old estimates and measurements from the primary sensors. If you operate the Kalman filter without primary sensors, it turns into a simple model-based extrapolator. Needless to say, extrapolating from the last known GPS position using steering and odometer data is much more accurate than a simple constant-speed extrapolation.
The correction stage of the Kalman filter becomes active whenever a measurement from a secondary sensor arrives. Each correction handles the measurement from only one sensor and updates the current estimates for position and orientation. There can be an arbitrary number of corrections with arbitrary or even irregular update rates (still no GPS in tunnels, therefore no corrections for several minutes).


Bandpass overshoots at the beginning and end of signal: why?

Well what you see is the transient response of your filter.  Your signal looks a bit a like a step signal at the beginning. At the end, you have a negative step, thus you will get another transient that will be the opposite of your first transient.

Plot the step response of your filter and compare the results.

Then ask yourself why does it respond that way?


16-APSK NDA Carrier Synchronization Problem

Symbol timing and carrier synchronization are two separate "tasks" in a demodulator.
There is lot of literrature on the internet to understand synchronization on a DVB-S2 signal.
But not sure that GNU Radio has the good ready block for it.
For symbol timing synchronization, there is Gardner timing error detector that suits well for QPSK modulation. An adapted version of it for 16APSK and 32APSK in a second order control loop could do the job.
For carrier synchronization, if you try to look at the 16APSK and 32APSK in polar representation (the constellation), the angle value modulo a constant gives you a phase error. Also in a second order control loop you can track the carrier this way. The modulo constant depends on the order of the modulation scheme you use. That is probably also why the one in GNU Radio may not work.
What I described before is called fine synchronization. Which means that it can be not enough to lock on large offsets.
You also need a coarse synchronization to help the fine one to lock. Depends on your requirements.
I hope this helps


It is important for the simulation that I'm running that there are only 4 possible modulation values before noise is added to the system.

So, you made a mistake when building your simulation, or your simulation already includes a frequency offset. The way you ask this, you probably did the former, but since we have zero knowledge of what exactly you've built: Can't help you with that, make sure your receiver simulation does exactly the same as what your transmitter simulation does. ("fix your program!")

I'll answer as if this was a system that's generally correctly designed:

  Thus, when using a QPSK modulator in GNU Radio I would expect to see, in the constellation diagram, a square with 4 points at each corner, but instead I see diamonds with thousands of points, even in the absence of any introduced noise.

That's not surprising! No two oscillators are exactly the same, and so your receiver and transmitter have a frequency offset.

A frequency is the linear change of phase over time. That means that your symbol phase at the receiver has a phase offset that grows constantly every time you look. Hence, the "diamond".

You need what you've been taught in your digital communications course: a frequency synchronization!

You can do that with a second order costas loop.


Speech Signals - Recognition vs Analysis

1) Questions of terminology like that are always debatable! I would say that the term "analysis" describes techniques to convert the speech signal into another representation, closer to the semantic level than a raw list of samples (such as a written transcription in the case of speech recognition, a visualization of the tone curve in the case of a system assisting Mandarin learner, the identity of the speaker in case of speaker identification system, a category of emotion in case of an emotion detection system) ; while "processing" describes problems in which both the input and the output are speech (speech denoising/enhancement, coding and decoding, source separation, 'creative' effects like pitch/timing modifications). Of course, there are contexts in which "processing" describes the whole field, analysis included.

2) Deterministic signal processing (filtering, transforms), statistical signal processing (auto-regressive models, sinusoidal model estimation), statistical modeling and machine learning, psychoacoustics,  knowledge of how humans produce speech... If you venture into recognition: language theory (finite state machines), natural language processing, some elements of distributed computing (for large model training), data structures/algorithms. If you venture into speech coding or enhancement applications: information theory, experiment design (for subjective evaluation campaign), DSPs/embedded systems with hardware acceleration or even in some cases hardware design. Tools and language: mostly matlab or scipy for research, mostly C/C++ for implementation; and a scripting language (perl, python) is always helpful when you have to deal with texts/annotations.


Using large DFT for OFDM to avoid the guard interval

Sure,  you can !
Practically, some systems are using OFDM system with number of subcarriers 4096, but the PAPR will be an issue in that case, and if your channel is fast-fading, the time-domain equalizer will be high-complex since one-tap frequency-domain equalizer will be useless.
Usually, choosing the length of the symbols depends on the channel delay, e.g., if your channel delay is 128,  you must use at least length of symbol 1024 and so on.


Matlab: What is the proper way to calculate mean square error for complex numbers

The mean square error is

e = ((abs(h_e1 - h_true(1))^2) + (abs(h_e2 - h_true(2))^2) + (abs(h_e3 - h_true(3))^2))/3;

But it is tediuos!

What if you had $1000$ terms instead of $3$? Do you want to add one by one?

Do it in the vector form for more clarity and ease of implementation:

$$\mathrm{mse}=\mathsf{E}(|\mathbf{h}-\hat{\mathbf{h}}|^2)$$

In MATLAB, you can use the function mse()

 h_true = [ 0.7071 + 0.7071i
            0.4243 + 0.4243i
            0.2121 + 0.2121i ];

 h_e = [ 0.8100 + 0.8100i
         0.5100 + 0.5100i
         0.1200 + 0.1200i ];

 e = mean(abs(h_e-h_true).^2); % or eqivalently e = mse(abs(h_e-h_true));


What happens when we reconstruct a signal after its Fourier transform is done and some frequencies are altered?

The spatial content is carried in the phase of the Fourier transform and not in magnitude. As Marcus pointed out, the difficulty with Fourier transform is that it doesn't allow us to immediately analyze the spatial information while still looking at the frequency. There are methods such as windowed FT which try to make that a bit better.

However, if you alter frequencies in the magnitude and reconstruct, you would get filtering effects. For instance, if we cut the low frequencies, than it would be a high-pass filter, and likewise, cutting high bands would cause a low-pass filter effect. This is analogous to applying, e.g. Mean/Gaussian filter to an image.

Regarding Wavelets, they try to address this time/frequency issue I mentioned above. Wavelets are limited in time and frequency. We slide the wavelets in time domain to construct the Wavelet domain. This provides resolution in the time domain.

On the application side, Wavelet denoising is shown to be superior to Fourier domain denoising. But, such superiority is application dependent, as for some applications one doesn't need that spatial analysis.

You could find some more information on the differentiation here.

Here is an example, based on this one, where we do a low-pass filter by keeping only the low frequency values (both positive and negative):

% make our noisy function
t = linspace(1,5,1024);
x = -(t-2).^2  + 2;
y = awgn(x,0.5);
F = fft(y,1024);

freqRange = 20; % range of frequencies we want to preserve

rectangle = zeros(size(F));
rectangle(1:freqRange+1) = 1;       % preserve low positive frequencies
y_half = ifft(F.*rectangle,1024);   % +ve low-pass filtered signal
rectangle(end-freqRange+1:end) = 1; % preserve low negative frequencies
y_rect = ifft(F.*rectangle,1024);   % full low-pass filtered signal

figure;
plot(t,y,'g--');
hold on, plot(t,x,'k','LineWidth',2);
hold on, plot(t,y_half,'b','LineWidth',2);
hold on, plot(t,y_rect,'r','LineWidth',2);
legend('noisy signal','true signal','positive low-pass','full low-pass','Location','southwest');

The resulting signals look like this:

You could see here that the altered frequencies effect the whole signal - wherever a high frequency existed gets smoothed. The spatial content is preserved in reconstructing the signal, because FT is a loss-less transformation, but we didn't get the chance to modify the exact spatial locations - the whole signal is effected by a change to the FT magnitude.


Implementing Convolution in Frequency Domain?

IDFT the smaller or both of the DFTs if needed.  Zero pad one or both of the kernel and image to make them the same dimension and size.  Re-DFT as needed, and now you can complex multiply the 2 DFT arrays element-by-element because they will now be the same size.


How do we distinguish between Image Recognition,Detection and Segmentation?

It can be a little confusing at times, and the terms are not completely independent.

Detection:

In detection, you are simply detection the presence of something. For example, you might design an algorithm to detect beach-balls from pictures of a beach. You would feed an image into your algorithm, and it will spit out an answer 'yes', if there was a beach ball, and 'no' if there was no beach ball. This is very simply what it means to 'detect' something.

Segmentation:

Segmentation is a method - a technique - an algorithm. So for example, if someone says "How did you design your beach-ball detection algorithm?", you might say "I used a segmentation algorithm". Or you might say, "I used a template-matching algorithm". Or you might say "I used a scale and shift-invariant matching score algorithm", etc. Those are all methods of doing something, and segmentation is just a method.

In particular, segmentation simply means a thresholding of your image into various areas. In a simple case, say you want to separate background from foreground. Then, you segment the image into two regions: A region that corresponds to the background, and a region that corresponds to the foreground.

Image Recognition:

Lastly, there is image recognition, which is simply, a very broad goal.  For example, you might say, "I achieved image recognition for beach scenes, by detecting beach balls using the method of segmentation, detecting water using the method of k-means, and detecting bikinis using the method of template matching.

Of course, those are just examples, but broadly speaking, this is how you should parse the different words and terminologies you come across in this field. Hope that helped!


I have finally found a great overview video on the topic from Richard Szeliski: "Reflections of Image-Based Modelling and Rendering"

The research started with stereo and depth estimation. The depth map can be used to "warp" image in 3D to give it some perceived depth. This does not create parallax effects though.

Another approach called Layered Depth Images used segmentation to several depth layers each having its own depth map. This allows parallax effect to occur, but extra care must be taken to matte foreground and background layers during rendering.

Fast 3D reconstruction of the sorrounding is possible using primitives (boxes, cylinders) matched to scene points. The 3D reconstruction is tedious and requires many images and extensive walk around.

Then, more image-based (rather than geometry-based) approach are presented, which are Lumigraphs and Light Fields. These allow some level of reconstruction of the plenoptic function and solve for rays crossing given viewing plane.

Finally, Concentric Mosaics is the no. 1 solution I have been looking for. A camera travelling around circle (tangential direction) is able to capture 3D plenoptic function of the scene. Each point thus can be viewed from different angles and the viewer can move inside the circle, see parallax, reflections, transparency and occlussion.

Extension is possible by adding support for rays in vertical direction - this will lead to capturing 4D plenoptic function from a camera moving on sphere (tangetial direction) or by using a spherical rig of many cameras.

Further advances in the field are expected with the introduction of plenoptic cameras (e.g. Lytro) that can capture the 4D plenoptic function directly. Some work on making plenoptic mosaics is already there (e.g. "Panorama Light Field Imaging").


Transformed Direct Form II filter: Z transform: reduce multiplications

It's instructive to look at the actual difference equation implemented by the system. If you define a sequence $w[n]$ as the output of the first (left-hand side) adder you get

$$w[n] = x[n] + Aw[n-1]\\
y[n] = Bw[n] + Bw[n-1]$$

Applying the $\mathcal{Z}$-transform gives

$$W(z) = X(z) + A \ W(z)z^{-1}\\
Y(z) = B \ W(z)\left(1 + z^{-1}\right)$$

which results in

$$W(z) = \frac{X(z)}{1-Az^{-1}}\\
Y(z) = \frac{B(1+z^{-1})}{1-Az^{-1}}X(z)$$

In the time-domain this last equation is equivalent to

$$y[n] = Ay[n-1]+B \left( x[n] + x[n-1] \right) \tag{1}$$

Equation (1) shows that this system can obviously be implemented with only two multiplications. If you want to use the Transposed Direct-Form II structure then just use the structure shown here (bottom figure), set $a_1=-A$, remove the branches with multipliers $a_2$ and $b_2$ and the lower delay element (because you only have a first order filter), remove the multipliers $b_0$ and $b_1$, and multiply the input signal directly with $B$.

EDIT: (in reaction to your updated question)

There are a few mistakes in your figure. Note that in the figure I linked to, the input is on the right-hand side. This might have caused some confusion. Anyway, in your figure where you multiply by $B$, just multiply by $A$ instead. Then remove the multiplication by $-A$ and add a multiplication by $B$ to the input signal $x[n]$, before the path splits. The node you refer to is simply an addition; they obviously forgot the 'plus' sign.

As for the equations, what I wanted to do is to get a difference equation for the output signal $y[n]$ in terms of $x[n]$ and of delayed versions of $x[n]$ and $y[n]$. This shows what the system actually does and in this way the description of the system becomes independent of the chosen filter structure. In order to do that, I used the $\mathcal{Z}$-transform which transforms difference equations into algebraic equations. So in the $\mathcal{Z}$-domain it is straightforward to eliminate the signal $w[n]$ and write $Y(z)$ directly in terms of $X(z)$. Transforming back to the time-domain yields the desired difference equation, which can then be implemented by any structure you like.


Concepts Surrounding First Two Samples in Difference Equation

This is just the transient effect of the filter due to zero initial conditions. The "solution" is to throw away the first two zero samples of $u[n]$ because they are not related to the input noise signal $v[n]$.


Signal that destabilises a positive feedback path with nonunity feedback

I am going to provide some initial information as simply as I can that should be helpful in the OP's attempt to get a basic understanding of control systems. Further study in Laplace transforms (and z-transforms for discrete-time systems) will be helpful in these frequency domain modelling approaches to control systems.
First understand the concept of a transfer function as the frequency response for a basic two port (one input / one output) linear system; the transfer function defines the magnitude and phase versus frequency of the output compared to the input.  This is often given in terms of the variable "s", but to avoid going down that path I am going to present it as a function of frequency only ($H(f)$).  Below shows an example system with a plot of the magnitude and phase versus frequency as given by the $H(f)$ equation. This shows us for any particular frequency at the input, what the amplitude of it will be at the output and what the phase shift will be between the input and output: every frequency will be transformed by the system by changing its amplitude and phase shift only, as long as it is a linear system no other frequencies will be created so this is a nice compact description that we will ultimately use to determine stability in a control system.

Now consider the simplest case of a system with feedback as diagrammed below with input given as $X$ and output given as $Y$ with the output fed back and subtracted from the input:

The relationship between $Y$ and $X$ is easy to show with simple algebra, with our goal to describe a new transfer function of this new closed-loop system. The new transfer function would be given as $Y/X$ (output divided by input). I don't show it below to not clutter the simplicity of the equation but $X$ and $Y$ are also functions of frequency as $X(f)$ and $Y(f)$, so would be the Fourier Transforms of time domain signals:
$$(X-Y)H(f) = Y$$
$$\frac{Y}{X} = \frac{H(f)}{1+H(f)}$$
The above equation would adequately describe our new closed-loop system as long as it follow the Nyquist Stability Criterion. To really show this completely to cover all cases, the OP needs to first understand the Laplace Transform and the generalized transfer function as a function of s (as in $H(s)$), but here in this simple introduction we can reduce it to a condition that is covered in most cases: that is if the gain of the system given by $|H(f)|$ is less than one before the phase of the system given by $\angle (H(f))$ exceeds 180°, the system will be stable and will have its own closed loop transfer function given by the expression above. We can see this through an iteration of a simple case with H(f) = -.8 at a particular frequency (gain = 0.8 and phase = 180°).
Consider the feedback system as an iteration with the input starting at 1: the output will be -.8, which fed back to the input as a subtraction will add 0.8 to the input, resulting in 1.8 at the output. This is then multiplied by -.8 resulting in 1.44 being added to the input resulting in 2.44 at the output...
With closer inspection we see that this is a geometric series:
$$1 + \alpha + \alpha^2 + \ldots = \sum_{n=0}^\infty(\alpha^n) = \frac{1}{1-\alpha}$$
Which converges only if $-1 < \alpha < 1$.
Where here $\alpha$ was 0.8 (due to the negative feedback in the loop) which would converge to $1/(1-.8) = 5.
If the gain of $H(f)$ exceeded one with the phase 180°, we would see how the same geometric series would not converge, as an example of an unstable loop.
This gives a very basic simplified explanation to help aid in an intuitive understanding, but further proper study would be to understand the Laplace Transform and complex systems, and the meaning of poles and zeros in the s-plane, and how a stable system must have all of its poles in the left-half plane, and ultimately the full details of the Nyquist Stability Criterion.


How harmonics happens when you mix two signals?

Amplitude modulation is a linear operation (excluding the carrier in the classical AM) and does not introduce harmonics due to the modulated signal. Sidebands are not harmonics. They are created as the message signal spectrum is shifted up to carrier frequency.
Frequency mixing of a message signal with a sinusoidal carrier, another definition of AM, does not introduce harmonics.
In a practical setting if there are some nonlinearities in the implementation (a deviation from the ideal case) you may have harmonics. Otherwise No.


What is the Fourier transform of $\delta(t-a)$

Your argumentation is perfectly right. The only thing you are missing is that the sum of complex exponentials becomes a comb function in the limit. See also this link and this matlab code:

N = 2048;
n = (0:(N-1))';

k = 4*(0:20);

signals = zeros(N,length(k));
for kk = 1:length(k)
    signals(:,kk) = exp(-1j*2*pi*k(kk).*n/N);
end

plot(abs(sum(signals, 2)));


To increase the power of x by 3 dB:

x .*= 10^(3/20)


How to transform lowpass fir filter to bandpass fir filter without using a built in function in matlab

Mathematically, applying a FIR with impulse response $h_\mathrm{lpf}[n]$ to digital signal is convolution:

$y = x * h_\mathrm{lpf}$, or thanks to the properties of the (discrete) Fourier transform,

$Y = X\cdot H_\mathrm{lpf}$, as convolution becomes multiplication.

Now, making a bandpass out of a low pass can be modeled by shifting the Frequency response $H_\mathrm{lpf}$ in frequency domain. "Shifting" can be represented by a convolution of the low pass filter with a dirac impulse at the desired center frequency:

$H_\mathrm{bpf}= H_\mathrm{lpf} * \delta_{f_\mathrm{center}} $

Again, convolution becomes multiplication when transformed to time domain. The inverse (discrete) Fourier transform of a dirac at ${f_\mathrm{center}}$ is a complex oscillation $e^{2\pi{f_\mathrm{center}}n}$, so this becomes

$y[n] = x * (h_\mathrm{lpf} \cdot e^{2\pi{f_\mathrm{center}}n})$


From the power spectrums you have plotted, it seems like there are a lot more deviations with the data. Perhaps when the song was recorded or encoded, the sampling rate (number of samples in each second) was less than the original one.


It is the error with respect to size of jump. It has nothing to do with error at discontinuity. As per the explanation given, suppose the perfect square is having values $0$ and $1$. Just after the discontinuity of a rise time, the perfect square would take value $1$, while the reconstructed wave would take the value $1.09$ at the peak of overshoot. Here the size of discontinuity is difference between perfect low value $0$ and perfect high value $1$ of the square wave. In general if it were $A$ and $B$, the overshoot percentage would have been
$$
\frac{(Ovs_{peak} - (B-A))}{(B-A)} \times 100
$$

The error of $0.5$ in the error graph is not related to overshoot. For perfect square wave of $0$ and $1$, at discontinuity it is either $0$ or $1$, but the reconstructed wave will be midway $0.5$ ($|0.5−0|$ or $|0.5−1|$). So max error will be at discontinuity and its absolute value will be $0.5$.


multiplications per second for FIR filter

To answer your questions:

The number of multiplies is the number of taps times the number of sample per second.  Given a sampling rate of 300MHz and 31 taps, you will have to do 9300 million multiplies per second.
As Paul R says, whether a multiply can be done in one cycle depends on the processor.  Some can, some can't, and some appear to by having a long pipeline - once the pipeline fills, the result of one multiply comes out per cycle.  You would have to find out what processor you are using to determine how many cycles it takes for each multiply.  Also note, this varies also with the type of number you are working with.  Integer multiplication is usually faster, but can cause you other difficulties (it is harder to do signal processing with integers.)
As Paul R says, really only the multiplies are of interest since most processors have a combined multiply and add instruction.  If you are working with a processor without that type of instruction, though, you will have to remember to count the adds when budgeting your clock cycles.  Since there are as many adds as multiplies when doing an FIR, knowing how many multiplies there are also tells you the number of adds.


Comparing a clean, noisy and enhanced signal

The answer to you question depends on how quality is defined in your specific application. Since you're working with audio signals, it is often beneficial to use the perceived quality to analyze the efficiency of an audio enhancement algorithm.
How to Determine the Perceived Quality
Methods to determine the perceived quality of an audio signal can be categorized into

subjective and
objective

methods. Subjective methods are based on human listeners that rate the audio signal according to some defined criteria. Objective methods use algorithms that use some measure, modeling parts of the human auditory system, to estimate the perceived quality.
Listening Tests
Usually, the most accurate way to determine the perceived quality is to perform a listening test and ask human listeners to rate the quality of different signals. Listening test methods comprise, e.g., the MUltiple Stimuli with Hidden Reference and Anchor and pairwise comparison approaches.
Objective Measures
A simple objective measure surely is the signal-to-noise ratio (SNR). In many cases, though, the SNR is only related to the actual perceived quality to a certain extent. Therefore, more involved objective measures have been created. Well-known ones are

the Perceptual Evaluation of Speech Quality (for speech signals) and
Perceptual Evaluation of Audio Quality (for audio/music signals)

algorithms.


PSD and $\lim_{T\rightarrow \infty} \frac 1 {2T} \int_{-T}^T x(t)\bar y(t)\,dt$

For power signals $x(t)$ and $y(t)$, the function

$$R_{xy}(\tau)=\lim_{T\rightarrow\infty}\frac{1}{2T}\int_{-T}^{T}x(t)\bar{y}(t+\tau)dt\tag{1}$$

is the cross-correlation of $x(t)$ and $y(t)$. So the expression you're asking about is the cross-correlation of $x(t)$ and $y(t)$ evaluated at lag $\tau=0$:

$$R_{xy}(0)=\lim_{T\rightarrow\infty}\frac{1}{2T}\int_{-T}^{T}x(t)\bar{y}(t)dt\tag{2}$$

Since the Fourier transform of $(1)$ is the cross-power spectrum $S_{xy}(\omega)$ of $x(t)$ and $y(t)$, the expression in $(2)$ is related to $S_{xy}(\omega)$ by:

$$R_{xy}(0)=\frac{1}{2\pi}\int_{-\infty}^{\infty}S_{xy}(\omega)d\omega\tag{3}$$

So the expression in your question is proportional to the integral of the cross-power spectrum of the power signals $x(t)$ and $y(t)$.


Are IQ-samples independent from the carrier frequency?

I'm not sure I understand your question, but here are some thoughts.

Say you have a passband (that is, with a carrier) quadrature signal: $$s(t) = x(t)\cos(2\pi f_c t) + y(t)\sin(2\pi f_c t).$$ Here $x(t)$ and $y(t)$ have bandwidth $B<<f_c$. If you sample this signal, the specific samples you get will obviously depend on $f_c$.

One almost never does that; for many reasons, it's better to remove the carrier. So, you can obtain the low-pass equivalent (complex envelope) signal $$s_{BB}(t) = x(t) + j y(t),$$ and sample that instead. Then, the samples will not depend on $f_c$.


Best Metric to Compare Sparsity of Vectors

I am sorry I cannot comment your answer due to my low reputation. Gini and your suggested sparsity ratio ($l_1(x)/l_2(x)$) both give me the same value for $\lambda$. But

The problem I still see is that I cannot take into account how well the vector is solving the equation $Ax-y$. I would like to combine the residuum $l_1(A\hat{x}-y)$ and the sparsity $l_1(\hat{x})/l_2(\hat{x})$ into a new metric value. Do you know something from the literature?
Gini only works with positive values, but my measurement vector has lots of negative values. In my eyes, this is a heavy disadvantage because many sparse vectors can have negative values.


How can I break a filter down into second order sections?

Below the results as sos matrix and gain. Every row in the matrix is a single second order section in the coefficient order of b0, b1, b2, a0, a1, a2. Obviously we need a0=1, and the sections are normalized so that b0=1 as well. The overall gain accumulated is in the single number "gain". The fact that this gain is actually in the order of 10^-9 shows that this is a numerically tricky beast.

Here is how this works.

Design the filter in poles and zeros (zp or zpk), NOT in polynomial form (ba). I imagine that the cheby1() function has an output option for that. If not, you can't use it.
The poles come in conjugate complex pairs, the zeros in this case are real, but you can think of them as a conjugate complex pair with a phase of 0.
Group one pole pair each with one zero pair. Start with the pole pair that has the highest magnitude and find that zero pair that's closest to it. This step is really important!!
Turn the poles and zeros into polynomial in 1/z of the shape (1-p/z)*(1-p'/z). Let's say the complex zero/pole is x+j*y. Then b0/a0 coefficients are always 1, the b1/a1 coefficients are -2*x and the b2/a2 coefficients are x^2+y^2.
Repeat over all poles. Each conjugate pole pair results in one second order section.

A word of caution, if I may: If none of things that I have listed above mean anything to you, than you probably should not be doing this. This is not trivial and many things can go wrong in non-obvious ways. Depending on what you are planning to do with the results, you really want to make sure that your code is fully vetted and properly tested. Unless you have a good idea on how to test this, I would not base anything mission critical on the results.

    sos =
   1.000000000000000   2.000000000000000   1.000000000000000   1.000000000000000  -1.824520683146258   0.835012167790563
   1.000000000000000   2.000000000000000   1.000000000000000   1.000000000000000  -1.825351261746257   0.854128332434601
   1.000000000000000   2.000000000000000   1.000000000000000   1.000000000000000  -1.845749373866514   0.899792317177862
   1.000000000000000   2.000000000000000   1.000000000000000   1.000000000000000  -1.890243419701735   0.963689345198141
   1.000000000000000  -2.000000000000000   1.000000000000000   1.000000000000000  -1.989176997083584   0.989218209927750
   1.000000000000000  -2.000000000000000   1.000000000000000   1.000000000000000  -1.996496707908275   0.996511870590091
   1.000000000000000  -2.000000000000000   1.000000000000000   1.000000000000000  -1.998714416931744   0.998722644911976
   1.000000000000000  -2.000000000000000   1.000000000000000   1.000000000000000  -1.999653030269727   0.999659262704001
gain =
    4.556034632804534e-009


what are the advanatages of using a scrambler

Scramblers are indeed used to "randomize" the data, i.e. to "whiten" the input sequence. The reason why this is necessary is that a transmission system has no control over the data the user is going to transmit. This causes problems because it violates the assumptions that are usually made when designing transmission systems, such as having independent data symbols. Note that usually the data are not random, as you seem to believe. A frequently occurring problem are long strings of zeros in the data, which can cause difficulties in timing recovery and adaptive equalization. These problematic sequences are removed (or, actually, made much less likely) by the scrambler.


I am also not an expert, but I would consider this: The 50 µm watermark is burned in the image and it is complete black (0 0 0) in rgb. It has a significant influence to any analysis oyu make on the image regarding statistics like mean, median, edge detection etc.

So it needs to be removed. The easiest method to remove it would be to find out if you microscope has a setting to metadata this information instead of burning it in the image.
If that is not available and you only have this image provided, then I personally would go as follows:

Segment the black "Noise" watermark area and
use that selection as a mask to fill with the mean color.
Do your analysis. Be critical of whether it influences the results and refine if necessary. (E.g. if it still interferes with edge detection add some random noise).

Either that, or, if you can afford it - exclude the area of the watermark altogether from the analysis.

But please correct me if anyone knows better.


One standard way to evaluate the quality of such techniques is to look at the distribution of correct clusterings versus incorrect clusterings.

This can be quantified by looking at the precision, recall and quality of the clusterings... but assumes that you have a "ground truth" (i.e. that you know the real cluster the data point belongs to).

Precision: $p = \frac{TP}{TP + FP}$, the percentage of positives that are true positives (and not mis-classified negatives).

Recall: $r = \frac{TP}{TP + FN}$, the percentage of true positives that were correctly identified.

Quality (AKA accuracy) = $q = \frac{TP + TN}{TP + TN + FP + FN}$, the percentage of correctly classified items from all classified items.

Here $TP$ is the number of true positives, $TN$ the number of true negatives, $FP$ the number of false positives, $FN$ is the number of false negatives.

If you have many clusters and need to evaluate their interaction, then a confusion matrix may be more useful.

If you want to go into this in even more detail, I'd suggest reading the paper:

  Andrew P. Bradley, The use of the area under the ROC curve in the evaluation of machine learning algorithms, Pattern Recognition, Volume 30, Issue 7, July 1997, Pages 1145-1159, ISSN 0031-3203, http://dx.doi.org/10.1016/S0031-3203(96)00142-2.
  (http://www.sciencedirect.com/science/article/pii/S0031320396001422)
  Keywords: The ROC curve; The area under the ROC curve (AUC); Accuracy measures; Cross-validation; Wilcoxon statistic; Standard error

A pre-publication copy of which is available.


Recursive form representation of FIR filter

The frequency sampling representation of a length $N$ FIR filter's transfer function

$$H(z)=\frac{1-z^{-N}}{N}\sum_{k=0}^{N-1}\frac{\tilde{H}[k]}{1-e^{j2\pi k/N}z^{-1}}\tag{1}$$

is a systematic way of finding a recursive implementation for any FIR filter. The coefficients $\tilde{H}[k]$ in $(1)$ are the DFT coefficients of the impulse response $h[n]$. I do not know of any other such systematic way.

However, there are special cases for which a recursive implementation can be found by using standard finite sum identities. One such class of finite impulse response transfer functions is

$$H(z)=K\sum_{n=0}^{N-1}a^nz^{-n}=K\frac{1-a^Nz^{-N}}{1-az^{-1}}\tag{2}$$

A special case of $(2)$ is the moving average filter

$$H(z)=\frac{1}{N}\sum_{n=0}^{N-1}z^{-n}=\frac{1}{N}\frac{1-z^{-N}}{1-z^{-1}}\tag{3}$$


Windowing before viewing a filter frequency response?

No, if you know the filter coefficients it's never a good idea to use a window. Why should it be? Instead of plotting the actual frequency response of the filter, you would plot the convolution of the filter's frequency response with the frequency response of the window, which would result in a "smeared out" frequency response. Especially for sharp frequency-selective filters the result would be misleading.

Note that the frequency response of an IIR filter is not computed by applying an FFT to its truncated impulse response (because that would be windowing with a rectangular window), but by computing the ratio of the FFTs of the numerator and the denominator. That's what Matlab's freqz does.

However, if you don't have the filter coefficients of an IIR filter, but just its impulse response, then you don't have much choice, and you have to use the impulse response. You must include as many terms as can be distinguished from noise and use an FFT. This corresponds to windowing with a (very long) rectangular window. I don't see much advantage of using any other window because you have to include practically all terms you can distinguish from noise in order to keep the truncation error small. See also the discussion in the comments below.


You can use corrcoef to compute the Correlation coefficients as a measure of linear dependence of the random vectors (i.e you transpose your row vectors, the variables are column vectors, and the rows are observations). With the vectors
u = [1, 2, 3, 4];
v1 = [10, 10, 10, 10];
v2 = [0, 1, 2, 3];

You have
>> R = corrcoef([u', v1', v2'])

R =

     1   NaN     1
   NaN   NaN   NaN
     1   NaN     1

There is a perfect match between u and v2 as v2' = u' - 1 showing that the two variables u and v2 are directly correlated. In the results above R is
$$
\mathbf R = \begin{pmatrix}
\frac{\operatorname{cov}(\mathbf u, \mathbf u)}{\sigma_{\mathbf u}^2}& \frac{\operatorname{cov}(\mathbf u, \mathbf v_1)}{\sigma_{\mathbf u}\sigma_{\mathbf v_1}} & \frac{\operatorname{cov}(\mathbf u, \mathbf v_2)}{\sigma_{\mathbf u}\sigma_{\mathbf v_2}}\\
\frac{\operatorname{cov}(\mathbf v_1, \mathbf u)}{\sigma_{\mathbf v_1}\sigma_{\mathbf u}}& \frac{\operatorname{cov}(\mathbf v_1, \mathbf v_1)}{\sigma_{\mathbf v_1}^2}  & \frac{\operatorname{cov}(\mathbf v_1, \mathbf v_2)}{\sigma_{\mathbf v_1}\sigma_{\mathbf v_2}}\\
\frac{\operatorname{cov}(\mathbf v_2, \mathbf u)}{\sigma_{\mathbf v_2}\sigma_{\mathbf u}}& \frac{\operatorname{cov}(\mathbf v_2, \mathbf v_1)}{\sigma_{\mathbf v_2}\sigma_{\mathbf v_1}} & \frac{\operatorname{cov}(\mathbf v_2, \mathbf v_2)}{\sigma_{\mathbf v_2}^2}
\end{pmatrix}
$$
The NaN results between (v1', u1'), (u', v1'), (v1', v1'), (v1', v2') and (v2', v1') is because v1' is constant, see MATLAB documentation on corrcoef, more specifically because
$$\begin{cases}
\operatorname{cov}(\mathbf v_1, \mathbf u) = \operatorname{cov}(\mathbf u, \mathbf v_1) = 0\\
\operatorname{cov}(\mathbf v_1, \mathbf v_2) = \operatorname{cov}(\mathbf v_2, \mathbf v_1) = 0\\
\operatorname{cov}(\mathbf v_1, \mathbf v_1) = 0
\end{cases}\quad\text{and}\quad\sigma_{\mathbf v_1} = 0
$$
In the end, each element of the matrix $\mathbf R$ is the covariance of two vectors divided by the product of their standard deviations.
EDIT:
Now let's say v1 isn't a constant vector and you have the three vectors as
u = [1, 2, 3, 4];
v1 = [11, 12, 11, 12];
v2 = [0, 1, 2, 3];

You then get
R =

    1.0000    0.4472    1.0000
    0.4472    1.0000    0.4472
    1.0000    0.4472    1.0000

Note that in the case here $\sigma_{\mathbf v_1} \neq 0$,  $\operatorname{cov}(\mathbf v_1, \mathbf u) \neq 0$, and $\operatorname{cov}(\mathbf v_1, \mathbf v_2) \neq 0$


What does photon noise mean in digital image processing or computer vision?

Ultimately at a fundamental level, an image signal as captured by a digital sensor is obtained by integrating the energy transmitted by the photons of the light wave.

If you are imaging a bright scene, there will be a large enough number of photons hitting the detector that this effect will not be detectable, but if the sensor is operating in a dark environment, individual photons may be detected and counted, and the number of photons that have hit a cell of the detector as a function of time is modeled as a Poisson random process.


Bandpass Filter window length: filter "in phase" with sinusoidal input

I find what you are saying really hard to understand. Let me try to decode it. Let me know if my decoding is correct.

  Say you have a bandpass filter and you place it on 40 discretely sampled points from a sinusoid with a period of 40.

This makes little sense.  What I think you mean is: I have a filter with impulse response $h(n)$ that I am going to apply to a sinusoidal signal of length 40 samples.

  Why is the bandpass filter (green curve on the image) only "in phase" with the sinusoid if you use a window length (impulse response length) for the filter of 1/4th the period of the sinusoid (here 40/4=10)?

Again, I'm finding it hard to fathom what you mean, but I interpret this as: If my impulse response $h(t)$ is of length 10 samples, then the output of my filter is a sinusoid that is in phase with the input.  When I use an impulse response of a different length (9 or 11) then the phase is incorrect.

For a start, your filter output will not generally have precisely the same phase as the filter input.  If your filter is a linear phase FIR filter, then as @MattL says, the filter will hvae a group delay of $(N-1)/2$ samples --- 4.5 samples in the case of a  length 10 filter.

Trying to reproduce your diagram (scilab code below) shows this effect: the green plot is the original, the red plot is the output of a length 9 filter (which has a delay of 4 samples), the black plot is the output of a length 10 filter (delay 4.5), and the blue is the output for a length 11 filter (delay of 5).

CODE BELOW

P = 40;
x = sin(2*%pi*[0:79]/P)

h10 = ones(1,10)/10;
h09 = ones(1,9)/9;
h11 = ones(1,11)/11;

y09 = filter(h09,1,x);
y10 = filter(h10,1,x);
y11 = filter(h11,1,x);

clf;
plot(x,'g');
plot(y10,'k');
plot(y09,'r');
plot(y11,'b');


This is very implementation specific and there is not one rule that one should be done before the other, and joint estimation another option where there isn't a clear order of events.

A consideration is how well your carrier tracking approach works under timing offsets, and similarly how well your code tracking approach works under carrier offsets. I recommend analyzing/simulating this specifically for your implementation since it is very application specific and you didn't provide any info about your waveform or algorithms for timing and carrier recovery. For example, the Gardner Timing Error Detector works very well with relatively large carrier offsets but has higher noise than the Mueller and Mueller detector (M&M), but the M&M being decision based does not perform well with relatively much smaller carrier offsets.

One suggested approach to quantifying this for your specific implementation/application: determine your pre-decision signal SNR with fixed SNR received signal input as you vary both carrier offset and timing offset. Also you can look at the noise of your symbol timing and carrier error estimators as another metric: determine what the noise is on the symbol timing error as your carrier offset is varied (When you are close to your track point for tracking considerations, and when you are further away for acquisition considerations). Similarly, determine what the noise is on the carrier error as your timing offset is varied. If you do this work, please do share it with us here.


Derivation of the LMMSE (Linear Minimum Mean Squared Error) Estimate and the MMSE Under Gaussian Prior

In the past I derived it as following:

It is a little different approach.
If it answers your question I will rewrite it in a proper LaTeX.
Regarding your question about the steps in the derivation you presented, it is using the Woodbury Matrix Identity (Both 12.26 to 12.27 and 12.28 to 12.29).
Related Answers:

Deriving the Matrix Inversion Lemma for RLS Equations vs the Woodbury Derivation.
Sequential Form of the Least Squares Estimator for Linear Least Squares Model.

In the answer above you may see some related derivations.


What is proof of Basic Hilbert Transform?

Hilbert transform has the definition that its CTFT is:

$$
H(\Omega) =
\begin{cases}
- \frac{\pi}{2}   &,& \text{ for } \Omega > 0 \\
+ \frac{\pi}{2}   &,& \text{ for } \Omega < 0 \\
\end{cases}
$$

from which you define the time-domian impulse response $h(t)$ associated Hilbert transformer from the inverse CTFT as:

$$ h(t) = \mathcal{IFT} \{ H(\Omega) \}  = \frac{1}{2\pi} \int_{-\infty}^{\infty} H(\Omega) e^{j \Omega t} d\Omega $$

You can evaluate the above integral or just use CTFT properties and pairs to conclude that:
$$ \boxed{ h(t) = \frac{1}{\pi t} } $$

Note that the associated impulse response $h(t)$ of the Hilbert transfomer is noncausal, two-sided and unstable and its Fourier transform can be expressed alternately as
$$ H(\Omega) = -j \text{sgn}(\Omega) $$
where $\text{sgn}(\cdot)$ is the sign function.

Hilbert transformer is most useful for defining and generating the analytic signal $x_+(t)$ associated with the signal $x(t)$, with has the property that

$$
X_+(\Omega) =
\begin{cases}
2 X(\Omega)  &,& \text{ for } \Omega > 0 \\
0            &,& \text{ for } \Omega < 0 \\
\end{cases}
$$


A multi-band filter is, as the name implies, a filter with multiple pass bands and stop bands. Here is an example of a multi-band filter design.


Using Kalman Filter to fuse two sensor readings of the same type

The answer to the first question is yes, it is possible to fuse data coming from the same type of sensor. For example, the fusion of two different accelerometers, each with a measurement quality, gives a weighted average where the weights are the inverse of individual variances (see Kays's[1] book for details).

The answer for the second is yes as well, it is possible to fuse IMU with GPS measurements. But note the one returns a position itself, while the other measures the second derivative of the position. The fusion demands well elaborated models. There are plenty of papers and books with models and its theory. I would recommend, for example, Farrel's[2] and Titterton's[3] books.You can easily find implementations of the fusion filters in internet. Open source projects of drones, for example, have it.

[1] KAY, Steven; Fundamentals of Statistical Signal Processing, Estimation Theory.

[2] FARREL, Jay A.; The Global Positioning System and Inertial Navigation.

[3] TITTERTON, D. H.; Strapdown Inertial Navigation Technology.


Your measurement of the symbol and noise energy is not correct.  Please note that the following is pseudo-code, not valid Python code.  EbN0_dB should be an input to this code, not an output.

Esymbol = sum(real(data*conj(data))) / numSymbolsInData
Eb = Es/numPayloadBitsPerSymbol
Eb_dB = 10*log10(Eb)
N0_dB = Eb_dB - EbN0_dB
N0 = (10.)**(N0_dB/10.)


How to quantify phase distortion of a filter specific to a signal?

For a phase distortion metric I recommend using “group delay variation”. The definition of Group Delay is the negative derivative of phase with respect to frequency. The Group Delay is the delay in time that “group” of signals over a band of frequencies would have. A frequency response that is linear in phase (constant group delay with no variation) is not considered a distortion since all the frequency components would have the same delay so no actual distortion of the signal results. When the phase is not linear versus frequency (as in your plots) different frequency components of the signal arrive at different times at the output of the system, which can result in considerable distortion. Group Delay variation can be quantified as peak variation, peak to peak, or rms as in any other distortion metric.

As MattL points out in the comments below, a more comprehensive metric would be deviation from linear phase, where linear phase strictly means proportional to frequency with no phase offset term. If a phase offset exists, it would not appear in the result for the group delay computation yet indeed contribute to a distortion due to a varying delay versus frequency (the Hilbert Transformer is an excellent example of this: in order for all frequency components to have a 90° relationship with the input waveform, each component must have a different delay). For further details on this see Matt's answer here:

Group Delay for Hilbert Transformer and Resulting Dispersion


How to get to shoulders localization on image photo?

Are there just shoulders in the video, or is there a head attached to them?

If it is the latter, then you can use vision.CascadeObjectDetector from the Computer Vision System Toolbox for Matlab.  It has the "upper body" mode, which effectively detects head and shoulders.


Please see the answer to this question How to implement Polyphase filter?
on how to implement a polyphase interpolator, and then note that the same structure without actually commutating the output can be used as a discrete variable delay; each output is the input within the passband of the filter at a different delay in uniform discrete steps. By stepping through the output continuously as we do in the interpolator, we achieve a higher sampled version of the input. By staying at one output we get the input at a specific delay based on that filter, by changing to another output we get another delay. The design of the filter as given in the link will properly create the filters such that each is a different delay of the other (with the same ampltitude response) in uniform steps.

Continuing with the example shown of "Resampling" further shows how a variable delay device results:

Here is shown a IS95 CDMA waveform with 2 samples per symbol (as indicated by the red dots on this eye diagram). The analog waveform that is represented by these samples is shown in blue.

If we used the polyphase structure as an interpolator, we would get all the samples in between the 2 samples per symbol we have, evenly spaced, given by the number of filter banks in the polyphase structure.

Therefore it should be clear that we can select any filter output that is closest to the desired delay we want to achieve (in this case within one symbol duration).

As mentioned in the comments a Farrow Filter is another approach to implementing a variable fractional delay in an FIR filter.


How can I decompose a signal into square waves?

What is described in the question is very near the Discrete Wavelet Transform (DWT) with the use of the Haar Wavelet.

The DWT decomposes a signal into a sum of dilated and translated orthogonal functions that do not necessarily have to be trigonometric. The DWT does not transform a signal from the time domain to a frequency domain but to a scale space where the "time" dimension is preserved. The Haar wavelet is effectively just one period of a square wave and due to its dilation and replication as the transform progresses it would appear as occurring at different frequencies. For more information on the link between decomposition level and frequency please see this link

Another transform that might be of help here, is the Walsh-Hadamard transform which does exactly that, decompose a signal into a sum of square waveforms which are orthogonal (Please note the sequency there as well).

For a brief example that seems to be near what you are after, please see this link

Hope this helps.


my Butterworth lowpass formulas do not agree with Fisher webpage

The problem is in the way you apply the bilinear transform. You have to use the appropriate (pre-)warping of the frequencies. Since the bilinear transform warps the frequency axis, you have to make sure that the corner frequency of the discrete-time filter is correct. One way to do that is as follows. The bilinear transform is defined as

$$s=k\frac{z-1}{z+1}\tag{1}$$

with some constant $k$ yet to be defined. If we denote the analog frequency by $\Omega$, and the normalized discrete-time frequency by $\omega$, Eq. (1) becomes for $s=j\Omega$ and $z=e^{j\omega}$

$$j\Omega=k\frac{e^{j\omega}-1}{e^{j\omega}+1}=k\frac{e^{j\omega/2}-e^{-j\omega/2}}{e^{j\omega/2}+e^{-j\omega/2}}=jk\tan(\omega/2)\tag{2}$$

Eq. (2) describes the frequency warping caused by the bilinear transform. If we use an analog lowpass filter with a normalized corner frequency $\Omega_0=1$, we must choose the constant $k$ such that for the desired discrete-time corner frequency $\omega_0$ the term on the right-hand side of Eq. (2) becomes $1$:

$$k=\frac{1}{\tan(\omega_0/2)}\tag{3}$$

where $\omega_0$ is the desired corner frequency of the digital filter. Eq. (3) and Eq. (1) define the appropriately normalized bilinear transform that you must use.

So for your example, the normalized first-order analog Butterworth lowpass transfer function is given by

$$H(s)=\frac{1}{1+s}\tag{4}$$

Applying the bilinear transform gives

$$H(z)=\frac{1}{1+k\frac{z-1}{z+1}}=\frac{z+1}{z(1+k)+1-k}=\frac{1}{1+k}\cdot\frac{1+z^{-1}}{1+\frac{1-k}{1+k}z^{-1}}\tag{5}$$

Ignoring the gain term $1/(1+k)$, the corresponding difference equation is

$$y[n]=x[n]+x[n-1]-\frac{1-k}{1+k}y[n-1]\tag{6}$$

With a desired corner frequency $\omega_0=2\pi/30$ you get from (3) $k=9.5144$ and $(1-k)/(1+k)=-0.80978$, just like on Fisher's website.


This approach with notch filter with not work. All clicks are impulse-like sounds and we know that an impulse has frequency content at almost every frequency. What you are trying to do, by applying the $6 \ \texttt{Hz}$ notch filter is to remove sinusoid with that particular frequency in your signal, isn't it? So do you hear that your clicks are such low-frequency sounds? Definitely not! Those clicks are spreading into whole of your spectrum (like very short burst of noise), and filtering out $6 \ \texttt{Hz}$ won't help!.
You also didn't shared any more informations about nature of your clicks and a signal (spectrogram, waveform or even recording would be great), because it might be the case that your signal is also impulse-like (i.e. recording of a shooting range), which will make the recognition and cleaning task bit more complicated. Also are your clicks the same over time? This could be also the case for digital systems, and simple correlation calculation should do the job. Additionally if you can sacrifice some of your recording quality and click's are very short, then you can use the median filtering. You can also use some time-frequency methods (based on wavelets, spectrogram calculation, LPC's) or  statistical approach, such as AR modelling, etc.
Here is some reading for you - I hope it will guide you and you will find appropriate method for your kind of signals.

Audio Click Removal using LPC
Godsil S., Rayner P. - Digital Audio Restoration - lot's of references
AR modelling based approach - book chapter
Detection of clicks using sinusoidal modelling for the confirmation of the clicks
Automatic detection and removal of impulsive noise in audio signals


Intuitively explain Bi-linear time frequency distributions, someone please?

As you mentioned - all these methods share common principle, they allow for representation of our signal in time-frequency domain. First thing to notice is that wavelets are very different from the STFT-like (linear) methods as they provide scale dependent frequency resolution (which in case of STFT is constant and governed by Heisenberg uncertainty principle). Probably you've seen this famous picture:

So called bi-linear transform mostly refers to Wigner Transform which can be considered as the Fourier Transform of signals' instantaneous autocorrelation function. It is governed by following equation:

$$W(t,f)=\int_{-\infty}^{\infty} x \left( t+ \frac{\tau}{2} \right) x^* \left( t-\frac{\tau}{2} \right) e^{-2\pi i f \tau} d\tau$$

where $^*$ denotes complex conjugate and $\tau$ is the time delay.

One might notice that in analogy to STFT, the shifted signal itself is considered as a window. It is also worth of mentioning, that Wigner Transform provides perfect time and frequency resolution at the same time - this is simply great for signals with fast changes of instantaneous frequency. On the other hand WT is the bi-linear transform, which implies major disadvantage: for each linear combination of two signals $x(t)$ and $y(t)$, result consists of two parts: auto-terms and cross-terms. This issue produces interferences in time-frequency representation known as "beats", which can make interpretation of WT difficult. If you are looking for more specific mathematical description, then please refer to article posted below - I see no point in rewriting equations.

Example is shown on a figure below, where two linear chirps were analysed. Cross-term is clearly visible in between, which makes interpretation of results troublesome. Sometimes it happens, that magnitude of cross-term is even higher than auto-terms.

Although some methods were described in literature how to deal with this problem. Please refer to the following article for more details.

  R.B. Pachori and P. Sircar - A new technique to reduce cross terms in the Wigner distribution


In filter design, why isn't it possible to have a frequency response with phase 0?

A zero-phase frequency response is real-valued, and, by a property of the Fourier transform, the corresponding impulse response is symmetric. Since it is symmetric (i.e. it is non-zero for $t<0$) it can't be causal, and consequently, it can't be implemented (without adding delay, which defeats the purpose).


Poisson summation formula and periodic summation of Fourier transforms

But we can choose the functions $x_1(t)$ and $x_2(t)$, even though they are equal at the sample points, in such a way that their transforms $X_1(f)$ and $X_2(f)$ are quite different from each other, and, therefore, so will their periodizations.

That's a wrong conclusion. If $x_1(nT)=x_2(nT)$ $\forall n$ then the discrete-time Fourier transforms of the respective sequences must be identical, and so must be their periodized spectra.

Note that $x_1(nT)=x_2(nT)$ and $x_1(t)\neq x_2(t)$ implies that at least one of the two functions is not sampled according to the sampling theorem, i.e., the function is not completely represented by its samples. Because otherwise equal samples would imply equal functions. This might be the source of your confusion. If both functions are band-limited and if the sampling theorem is satisfied then $x_1(nT)=x_2(nT)$ implies $x_1(t)=x_2(t)$.


find a line go through the center mass and that line must be inside object area

How about:

Separate the foreground from the background and keep a list of indices or a logical mask from the background
Convert the image to a (rho, theta) coordinate space, meaning that each pixel is encoded by the angle and length of vector from the origin to itself. You can scan along this vector to see if you encounter any background pixels. Have a look at the so-called Trace Transform
Alternatively you can stay in the original coordinate domain and use Bresenham's line algorithm to get all pixels belonging to a scan-line between a point and the center and then comparing them against the list of background labeled pixels.

Perhaps you can optimize a bit by starting on the edge of the foreground shape, working your way inward to the center, discarding any scan-line the moment it encounters a background labeled pixel.


What is the maximum frequency that can be captured by a given sampling rate?

The maximum frequency is also limited by how long you sample at a given sample rate.  If you only sample for a finite length of time at 8000 sps, the bandwidth will be less than 4kHz.  A 4000 Hz signal (or any integer multiple thereof) will be ambiguous, or simply disappear.

A signal somewhat below the Nyquist frequency (depending on how long you sample at a given S/N or accuracy) can be "captured", as well as any signals between two successive integer multiples of Fs/2 given a-priori knowledge of no aliasing from any other spectra inside any other folds.


Explanation of fundamental filtering's consequences on signal

As already mentioned by @Fat32 in the comments to your question what you're observing is not commonly referred to as overshoot.

As there are no units on both plots comparing the results in both figures is somewhat difficult. The following assumes that the horizontal axes in both plots are just the indices of either the time or frequency domain values (otherwise the graphs cannot be directly related without additional information and/or do not make sense).

The exact effects of what you're observing will generally depend on the frequency, amplitude and relative phase of the harmonics present in your signal.

A rough visual comparison of the green and blue lines in the first plot indicates that indeed the first harmonic has been removed. The difference between the green and blue lines is roughly a sine wave (phase zero) that fits into the window (i.e. no difference at the start, end and 1/2 of the window; largest difference at roughly 1/4 and 3/4 of the window).

Note that if the amplitude and/or phase of the first harmonic are different the time domain signal, resulting after subtraction of the first harmonic, will probably look completely different.


Do you know Calculus?

All those shapes correspond to bounded values over an interval of the function or its first and second derivatives.

There are many ways to approximate derivatives in discrete sequences.

You can design a metric for each of your conditions based on the derivatives, then evaluate those metrics on a sliding basis.

1) A step occurs when the first derivative is relatively close to zero, then has a spike, then returns to zero.

2) Same as 1, opposite sign.

3) An upward line exists when the first derivative is bound near a positive value and the second derivative is bound near zero.

4) Same as 3, opposite sign

5) A hill (approximated by a parabola) occurs when the second derivative is bound near a negative value.

6) Same as 5, opposite sign


Where is the "time and frequency" in a Discrete Wavelet Transform?

DWT is motivated by the ability to analyse the signal in time and frequency with an adaptative time and frequency resolution.

Some clever algorithm allows you to perform a wavelet analysis with a recursive implementation. Matlab's DWT does that, matlab's dwt perform only one level of DWT, therefore, A is only the output of the low-pass filter while D is the output of the high-pass filter. There is no time and frequency here.
It only gives you a very coarse information about the time-frequency content of the signal.
Please check how to implement a DWT.

Second question : To see the novelty, you have to compare that to short-time frequency transform. With the DWT the time-frequency resolution depend on the frequency while with the STFT it's always the same. This is a VERY DESIRABLE property : You want a short time windows for quicly mooving high frequency signals and a long window for slower low frequency signals.

Moreover, please consider that the filtering/downsampling implementation is a very refined and efficient one. DWT is an incredible tool for analysis, denoising, feature extraction etc ...


HINT:

If you have a sequence

$$\hat{x}[n]=\begin{cases}x[n],&n\text{ even}\\0,&n\text{ odd}\end{cases}$$

then the DTFT of $x[2n]$ can be written as

$$\sum_{n=-\infty}^{\infty}x[2n]e^{-jn\omega}=\sum_{n=-\infty}^{\infty}\hat{x}[n]e^{-jn\omega/2}=\hat{X}(e^{j\omega/2})\tag{1}$$

You can obtain $\hat{x}[n]$ as

$$\hat{x}[n]=\frac12\left(x[n]+(-1)^nx[n]\right)\tag{2}$$

Taking the DTFT of $(2)$ to express $\hat{X}(e^{j\omega})$ in terms of $X(e^{j\omega})$, and using $(1)$ will give the desired result.


Designing a discrete-time inverse system

Correct answer is in the comment already but for completion: Your transfer function is

$$H(z) = b_0 + b_1 \cdot z^{-1}, \mbox{ with } b_0 = 1, b_1 = .9$$
which makes the inverse
$$H^{-1}(z) = \frac{1}{b_0 + b_1 \cdot z^{-1}} $$
which corresponds to the difference equation
$$y(n) + 0.9 \cdot y(n-1) = x(n) $$ or
$$y(n) = x(n) -0.9 \cdot y(n-1)  $$


How do I deconstruct a signal with fast fourier transformation in Python?

If you have a periodic signal, that is a repeating waveform, you can most easily find the constituent component parameters by aligning the Discrete Fourier Transform (DFT) frame on a whole number of cycles.  Suppose you chose four.  Then bins 0, 4, 8, 12, etc. contain the relevant values.

For a signal:

$$ x[n] = A \cos( \omega n + \phi ) $$

If

$$ \omega = \frac{2\pi}{N} k $$

where $k$ (cycles per frame) is an integer, then the $k$th bin will have a value of

$$ X[k] = \frac{A}{2} e^{i\phi} $$

So, you can read the tone's parameters right from the magnitude and angle of the DFT bin.  (FFT is an implementation of DFT)

If you aren't aligned on a whole number of bins, there is an interative process you can do with the formulas found in my blog articles.  This takes considerably more calculations though.


It is easily possible to design FIR filters with a complex frequency response specification, i.e., with a desired frequency response

$$D(\omega)=|D(\omega)|e^{j\phi(\omega)}\tag{1}$$

where both the desired magnitude $|D(\omega)|$ as well as the desired phase $\phi(\omega)$ can be prescribed. A least squares approximation of $(1)$ by an FIR filter just involves the solution of a system of linear equations. I've written a Matlab/Octave function solving that problem: lslevin.m.

Note that it is important to choose the specification $(1)$ in such a way that it can be realistically achieved by an FIR filter of a given length $N$. E.g., the average delay implied by the desired phase $\phi(\omega)$ should be close to $(N-1)/2$, which is the delay of a linear phase FIR filter of length $N$.


Oversampling while maintaining noise PSD

One important point you may be missing is assuming that the noise of the ADC is truly white over your sampling bandwidth- it isn't! You will have spurs from signal correlations and ADC non-linearities that will limit the noise floor you can achieve. This will be specified by the SFDR (spurious-free-dynamic-range) of the ADC for guidance as to what you may be able to achieve. Since the location of the spurs will be difficult to predict, you should assume that the strongest spur could be within your bandwidth of interest, in which case no linear filtering will reduce the power of the spur which exists at a single point in frequency.

That said, the filter will reduce white noise components, and the best metric is the EQNBW (equivalent noise bandwidth) for your filter; which is the equivalent bandwidth of a brickwall filter that would remove the same amount of noise. In this case the SNR improvement (processing gain) is 10Log10(2EQNBW/fs) assuming a low pass filter where EQNBW is the bandwidth from DC to the equivalent cutoff frequency.


Regarding Bode plots; $H(s)$ and $H(j\omega)$

Note that each complex pole $s_{\infty}=\sigma + j\omega$ of the transfer function $H(s)$ contributes to the system's impulse response a complex exponential of the form

$$e^{s_{\infty}t}=e^{\sigma t}e^{j\omega t}$$

The term $e^{\sigma t}$ is real-valued and represents exponential damping (assuming that the system is causal and stable, i.e. $\sigma<0$), and the complex term $e^{j\omega t}$ represents an oscillation at frequency $\omega$. So $\sigma$ is a damping constant.


Purpose of multiplying fft output by Hann window?

Can anyone briefly describe the DSP theory behind multiplying a flat (rectangle) spectra with a Hann window of the same size?

Without any reference to the actual model and its purpose, the only thing one can say about this is that due to the convolution theorem, the described multiplication results in some form of a digital filter. That is all one can say about "modulating" the frequency domain with any sort of function.

BUT, based on this excellent exposition of the way parabolic equations can be solved via the Discrete Fourier Transform, it seems that the multiplication is performed in both the spatial and frequency domain to isolate specific parts of the signal "gracefully".

Now, the authors seem to indicate that the multiplication is done to "...avoid aliasing...", but I cannot see how the results of these operations can lead to aliasing when:

One is in full control of everything about the simulation
The first "plane" (the one that is generated at the source) would only contain harmonics form the desired spectrum and be zero elsewhere anyway.

So, I cannot see how extra harmonics are likely to creep in between integration steps to have a need to isolate them.

Therefore, a more likely reason might be to be able to "isolate" (or focus on) only a part of the domain in a graceful way, rather than setting it to zero.

Hope this helps.


Matlab: Separate data samples from signal

This is a little hacky, but it uses the three models you have to find an "eigen"-pulse and then takes successive inner products of a normalized signal segment with this "eigen"-pulse to give a measure between -1.0 and 1.0.

The red line in the image is your data.  The black line is this normalized measure.  You'll have to sanity check it (no detections within 2 samples of a previous detection ?), but other than that selecting anything above 0.8 should let you detect the start of the pulses.

R Code Below

  #31201

  arm <- c(12,-23,1,-2,-40,-45,80,85,14,-28,10,-22,-100,-105,80,85,14,-2, 16,-2,-200,-45,305,85,1,-2,1,-2)

  #1 [-40,-45,80,85]
  #2 [-100,-105,80,85]
  #3 [-200,-45,305,85]
  rep1 <- c(-40,-45,80,85)
  rep2 <- c(-100,-105,80,85)
  rep3 <- c(-200,-45,305,85)

  repMat <- rep1 %*% t(rep1) + rep2 %*% t(rep2) + rep3 %*% t(rep3);
  mean1 <- (rep1+rep2+rep2)/3

  eigs <- eigen(repMat)

  plot(rep1, ylim=c(-350, 350), col="red", lwd=6, type="l")
  lines(rep2, col="green", lwd=4)
  lines(rep3, col="blue", lwd=2)
  lines(mean1)
  lines(eigs$vectors[,1]*sqrt(eigs$values[1]), col="purple", lwd=1)

  eigs$vectors[,1] %*% rep1

  similarity <- rep(0,length(arm))
  for (idx in seq(1,length(arm)-3))
  {
    patch <- arm[idx + seq(0,3)]
    patch <- patch / sqrt(sum(patch ^ 2))
    similarity[idx] <- (eigs$vectors[,1] %*% patch)
  }

  plot(similarity, type="l")
  lines(arm/max(arm), col="red")


Scatter plot: calculate box where 80 % of the points are

Your question seems to imply that there will always be just a single cluster of points.  If this is the case, I would go with a fairly simple heuristic like the following-

Calculate the mean and variance of the point positions in both the x and y directions.  The mean is the center point of your box, and you can use the ratio of the x/y variances to determine in which directions the box should grow faster.  For instance, if the x variance is "1" and the y variance is "2", I would have the box expand twice as fast in the y direction.

Now simply expand the box until you have at least 80% of the points contained.

If there is ever more than one cluster of points you may need to use a more sophisticated solution that uses clustering algorithms to identify the clusters and a box for each cluster.


A good stringed instrument music pitch detection/estimation algorithm will do the opposite, e.g. it will not ignore overtones.  Instead it will pay attention to the harmonics, specifically the harmonic train and its spacing, as this is a stronger indication of human perceived pitch than spectral content at the fundamental frequency (which could be nearly or completely missing, depending on the microphone and channel characteristics).

If you want to work in the frequency domain, than look into cepstrum/cepstral analysis and/or the Harmonic Product Spectrum algorithm, which can estimate pitch from spectral periodicity.  Or you can instead work in the time domain using weighted/interpolated autocorrelation, ASDF, AMDF, etc. type algorithms for pitch estimation.

If you slice your time windows too small, then there will be very few full cycles at he fundamental pitch frequency, which causes FFT frequency estimation to have very poor resolution.  A parametric estimator within and across windows may work better than a bare FFT.


Confused in Smoothing a signal in the Frequency domain

1) Is this a theoretically accepted method to smooth data in DFT domain or only applies in time domain ?

Sort of. It certainly can be done this way. Doing it directly in the DFT (instead of the PSD) is risky since it's a complex number you can get a lot of cancellation if the phase is fluctuating a lot.

You also need to decide how you manage the "edges": do you do circular or linear convolution

  2) Are there different methods to smooth DFT data?

Yes. There are as many methods as there are different sets of requirements and application. You plotted your data on a log frequency scale, so you may be better off with "1/Nth octave smoothing algorithm"


Isosurfaces from three dimensional column data: methods

You are given a series of $N$ triplets of values ($P_1$, $P_2$, $P_3$), pertaining to physical measurements. $N$ is in the order of $100$ to several thousands. My colleague would  like to draw isosurfaces from this dataset.

This is a little bit broad but perhaps still with some chances of accurate answers.

To an extent, it is important to know where the data comes from because it affects how the iso-surfaces are reconstructed. A Digital Elevation Model, for example, is essentially a 2D construction but volumetric data are inherently 3D. In the latter case, a trully 3D surface has to be fitted through the points of the contour which is much more challenging.

In both cases, the marching cubes algorithm is the standard way to extract the iso-surface but other methods involving a shrinking elastic surface (basically, 3D generalisations of active contours) have also been used.

In terms of practically applying it:

For relatively small point-clouds (a few hundreds of points), Blender with Point Cloud Skinner, can produce quick results. I have used this to create a surface through Magnetoencephalography sensors (whose locations are known) which ends up looking like a helmet. If your point cloud describes a convex shape, this will have no problems recovering it. I think that it works similarly to QHull (i.e, extracts the convex hull) but it has a threshold parameter that to an extent allows it to account for holes.
Paraview. Incredibly useful and with Python scripting capabilities. The software handles loading the data and several standard operations and allows extracting the results. The good thing about Paraview is that it includes a "raw import" functionality which you can use to load data in whatever format you may have it in. Isosurfaces are produced with "Contour". This is straightforward for volumetric data. If you already have the points, there is a "Table to structural grid" function that will bring the data to the expected form prior to applying "Contours".

I would consider this the faster option over the Visualisation Toolkit but if you need more control over the result, you might find the VTK more useful.

Meshlab. Meshlab is specifically built for handling huge point clouds including fitting surfaces. Meshlab is interesting because there are various scripts you can use around it, depending on what you are trying to achieve.
The Point Cloud Library and its Python bindings. This is similar to Meshlab in that it can handle huge volumes of data and it includes a surface fitting capability. It is however more challenging in setting up and operating.

Hope this helps.


It symbolizes the input-output characteristic of a half-wave rectifier, which only passes positive input signals:

$$y(t)=\begin{cases}x(t),&x(t)>0\\0,&x(t)\le 0\end{cases}$$


Discrete convolution is a linear operation, so yes

$$\begin{align*}x[n]*y[n] &= (x*y)[n]\\
\\
&= \sum_{m=-\infty}^{\infty} x[m]y[n-m] \\
\\
&= \sum_{m=-\infty}^{\infty} x[m]\Re(y[n-m])+ jx[m]\Im(y[n-m]) \\
\\
&= \left[\sum_{m=-\infty}^{\infty} x[m]\Re(y[n-m])\right] + \left[\sum_{m=-\infty}^{\infty} jx[m]\Im(y[n-m])\right] \\
\\
&= x[n]*\Re(y[n])+jx[n]*\Im(y[n])\\
\end{align*}$$


Even more on: Kalman filter for position and velocity

There are some problems:

  1-Did I inject the noise at the correct location? Is it appropriate to
  inject it prior to the double integration?

The noise you injected accounts for process noise at the acceleration state, which you seem to assume as constant plus random disturbance. But then your true state generation is wrong.

  2-What's a good real-world noise power level to use in this
  application? I picked a number at random, but I imagine whatever I
  pick has to be considered in my Q or R matrices. I know if I bump the
  noise up to something like 50, the filter output and measurement value
  go wild, but still track each other.

Which noise? Measurement or Process noise? Where is your measurement noise added? On the simulink model or on the script file?

  3-To get the simulation to work I had to change my H matrix to 1, but
  [1 0] was specified in the example. Does anyone know why?

H =[1 0] is correct if your state $X = [x ~~\dot{x}]^T$ has 2 elements such as position $x$ and velocity $\dot{x}$. Becasue your measuring just the first state $x$, the position.


Well, it's important to understand that face recognition (just like object recognition process) is a two-stage process, the first one is the Computer Vision phase, which is to represent the image (face) in feature space that works well (depending on the task you want to do) so for example, the image itself is a set of ordered pixels (let's say 300x300), which makes each image representative by 90000 numbers, this is the dimensionality of your data points (images) and just FYI, 90000 is insanely huge, if you try to use Neural Networks, Support Vector Machines, or any other classification tool using this high dimensions, you will absolutely learn nothing - this is the curse of dimensionality (you will need to provide tremendously larger set for training, in other words, to learn any classifier you will need many more images than 90000, which is hard and unreasonable thing to do).

So you will need to represent the images in a different space (called feature space), for face recognition, there is something called eigenfaces (look it up), but the idea is that you find the most discriminate features of the faces then for each image, you convert that image to a set of numbers that will use the discriminant features to describe the image, so each image will be represented as a set of numbers (.3, .1, .1, ...) which will be much smaller than 90000. By the way (.3, .1, .1, ...) will mean something like for example .3 there's a difference between hair and face, .1 there's difference between face and background and so... However, this is a really hard step (feature selection) and it may vary based on the kind of images you have (like the basic eigenfaces only work when your images are aligned, are your images aligned?). Of course using Color Histograms is another approach to represent the image/face (personally (this is only my opinion), I think that color histograms won't be so efficient for face recognition, I might be wrong and would love to know if anyone knows)

The second step is the classification step, here you can use Neural Networks (which also needs a lot of parameter tuning), Support Vector Machines (popular), or decision trees (simple and interpretable), adaboost (popular), and some other approaches. - so you will provide your images in the feature space instead of the pixel-representative and hopefully you'll learn a good classifier that will classify with high accuracy.

Hope that helps


The two types of stability and "Why exponential"

Let's consider discrete-time systems. The generalization to continuous-time systems is straightforward.

First note that BIBO stability does not require that the impulse response be square summable. It is necessary and sufficient that the impulse response be absolutely summable:

$$\sum_{n=-\infty}^{\infty}\left|h[n]\right|<\infty\tag{1}$$

In the context of stability of a causal LTI system, an exponentially bounded impulse response must satisfy

$$|h[n]|\le a\rho^{n},\quad n\ge 0\tag{2}$$

for some positive constants $a$ and $\rho$, with $0<\rho<1$.

If a causal LTI system is described by a rational transfer function, it can be shown that condition (1) for BIBO stability is equivalent to the condition that all poles of the transfer function lie inside the unit circle of the $z$-plane. This is equivalent with the system's impulse response being a weighted sum of decaying exponentials, and, consequently, condition (2) is always met by choosing $\rho$ equal to the largest pole radius, and by selecting an appropriate value for $a$.

In sum, condition (2) is a special case of general BIBO stability for LTI systems with rational transfer functions.

Also take a look at this answer, which also mentions non-linear systems.


In your example, the sinusoids are still orthogonal to one another; the key is that you must integrate over the correct period when calculating $P$. All of the sinusoids have frequencies that are multiples of $\frac{\phi}{T}$, so their resulting sum will have period $\frac{T}{\phi}$.

You should therefore take this fact into account when calculating the average power over a period of $s(t)$:

$$
P = \frac{\phi}{T} \sum_{k = 0}^{N-1} \sum_{k' = 0}^{N-1} a_k a^*_{k'} \int_{0}^{\frac{T}{\phi}} e^{j2\pi \phi (k - k')t/T} dt
$$

With this formulation, then you do get the useful property that the integral is zero for all $k \ne k'$, which allows you to simplify the expression to be simply the sum of the powers of each of the sinusoids.


In OFDM, what is the length of the cyclic suffix that is added before windowing?

You, as the designer of an OFDM system, define the length of the CP based on your requirements.

The main requirements are:

Long enough to "swallow" the duration of the longest expected channel impulse response to inhibit ISI / make the linear convolution channel effect look like a circular convolution
Not longer than necessary, since it "wastes" transmit time, and thus reduces data rate
(optional, depending) in a length that allows for efficient synchronization (e.g. via Schmidl&Cox)

So, there's no single answer to this. The answer is:

  It depends, and when you've fully understood the reason why we use CP in OFDM systems, this will become very much more intuitive.


Fourier transform for non-sinusoidal signal

The Fourier Transform will decompose your non-sinusoidal signal into harmonics, dominantly odd harmonics since your distortion appears symmetrical, and the amplitude as you derive would be the amplitude of the relevant sinusoidal harmonic (so in your case it looks like the fundamental is shown, so we are seeing the amplitude of the first fundamental harmonic which is a sinusoidal signal). It will be less than the peak amplitude shown in the time domain plot, which is the composite of all the harmonics.

Other factors that will effect the amplitude is spectral leakage due to your use of a rectangular window (so to the extent the harmonics fall into a sidelobe of the kernel for your rectangular window) and scalloping loss (to the extent the fundamental frequency is between an integer sub-multiple of your sampling frequency. Both of these effects are described in more detail in my favorite paper by fred harris: fred harris On the Use of Windowing

I am surprised that we do not see any evidence of these harmonics in your frequency plot but this may be because you are only showing a portion of all the frequencies or perhaps because the magnitude is not on a log scale we are not able to make out the harmonics. For example, there should definitely be a signal at three times the fundamental shown. The fundamental appears to be close to 0.04; meaning we should see a harmonic at 0.12; if your frequency of 0.1 represents the Nyquist boundary, then the image of this would be at 0.08 assuming a real signal as you have plotted. (Could you please udpdate your plot to have a log (dB) magnitude scale so we can be sure nothing else is astray?).

That said, we do see that your amplitudes are actually quite close visually from the plot (approx 0.17 in the FFT vs approx 0.19 in the time domain plot), so you must be concerned with the small difference, which is accounted for from the effects described above (the 0.17 is the amplitude of the primary tone as modified by the spectral leakage of your rectangular window and any scalloping loss).


Why does the amplitude of a discrete convolution depend on the time step?

Your computer doesn't compute the continuous integral, it does discrete convolution, which is just a sum of products at each time step.  When you increase dt, you get more points in each signal vector, which increases the sum at each time step. You must normalize the result of conv() according to the length of the vectors involved.


The symbol duration $T_p$ is related to the symbol rate $R_s$ by the simple relationship $$T_p = \frac{1}{R_p}$$
In the case of BPSK, $R_p$ is equal to the bit rate $R_b$, and $$T_p = \frac{1}{R_b}$$
For constellations with $M$ elements, each symbol encodes $k = \log_2(M)$ bits, the symbol rate is $R_p = R_b / k$, and the symbol duration is $$T_p = \frac{1}{R_p} = \frac{k}{R_b}$$