The transient edge or pulse shape of a tone burst represents a modulation.  The resulting modulated signal will contain frequency components higher than the tone frequency.  Thus you may need a sample rate much higher than just over 2X the tone frequency to reproduce a waveform accurately enough for cross-correlation to match the edges of the pulse envelopes.  This is likely true for your low tones and for voices, where the sample rate can be 10X to 200X higher than the tone or forment frequencies.

Filtering out frequencies just below the tone carrier can also remove any lower sidebands of the modulation, thus likely increasing the error in envelope matching by cross-correlation.


why do low pass and high pass filters generate a phase shift?

Because they need a certain time delay in order to "decide" what the output should look like if the output changes a lot. Imagine you see a number series and your brain has to decide what the mean value of these number series is. You will need to look at several numbers to arrive at some reasonable conclusion (let's assume the numbers are represented by a noise distribution with the same mean). The same applies to filters that weigh signal values to arrive at the filtered output signal.

Any filter that does not employ prediction of output will have phase shift. The only way to get rid of phase shift completely (or keep it to fractions of dB) is to use a predictive filter - for example kalman filter. Such a filter, in the simplest form will estimate the derivative of the signal and use that to filter the original signal in a predictive fashion. If you have a model of your signal you can do even better. Since prediction and sensor measurements can be processed in the same time step, from a discrete time point of view you can have zero phase shift if prediction is able to estimate the original signal precisely.


I understand they  implement the convolution by Lo and Hi db2 filters given after the sentence:

  Rather, we propose to use db2 filters. The filter coefficients for the
  db2 filter are represented as

I have not seen in the paper a claim that Lo_D and Hi_D in Figures 5b and 5c should be half-size. They are low-pass and high-pass branches. However, they should be followed by a downsampling to implement a  one level-DWT, or at least something similar.

There are different ways for  implementing a one level of a standard critically sampled DWT: direct, lattice, lifting, polyphase, Fourier. They should end up all with the  same number of samples as the original signal (away from the borders).


How much can we downsample by? and how much can we upsample by?

Nyquist. If your signal is real and has a BW of B = 100 Hz (meaning DC to 100 Hz), then the most you can downsample by is 2, in order to obey Nyquist and have a sampling rate that is at least 2B = 200 Hz.  Note you will typically want to sample more beyond this limit to allow room for reasonable filtering.


Why are Gold codes and Kasami codes used instead of pure m-sequences?

There are very few m-sequences of any given length with good cross-correlation properties.  Their autocorrelation properties are excellent, but the cross-correlation properties are variable.  For example, there are 18 m-sequences of period 127 but to have good cross-correlation properties, one must choose a set of no more than $\require{cancel}\cancel{\text{3}}$ 6 m-sequences out of the 18 available. As a consequence, when several sequences with good cross-correlation properties are needed, one needs to use Gold sequence sets  or Kasami sequence sets which have good cross-correlation properties and good (though not excellent) autocorrelation properties.

For more than you probably want to know, I refer you to the paper D.V. Sarwate and M.B. Pursley, "Cross-correlation properties of pseudorandom and related sequences," Proc. IEEE, vol.68, pp.593-619, May 1980. It includes a detailed discussion of the Gold sequences as well as the small sets and large sets of Kasami sequences.


DFT - Removing window effect in spectral domain with convolution

There is, in fact, a downside to your suggestion: the signals you showed are all very clearly divided into their frequency components, but in general, real-life signals tend to be more noisy.

Depending on the application, you'd want as much leakage dampening (higher main lobes/smaller side lobes of a signal frequency in the transformed spectrum) or, as another example, the narrowest possible main lobe.

In your plots, it is visible that smoothing the magnitude spectra with a window does the opposite of these: the main lobes get wider and smaller, while leakage products from transforming a finite time signal gain in power. Applied to noisy signals, this would yield a significant disadvantage.

Your suggestion, however, is still pretty useful for peak identification!


How to effectively teach topic of frequency domain?

In my personal opinion, it's best to study how the master teachers of the field have presented this information, and to the extent possible, copy their approaches and their style.

There are two teachers in the field who have more or less taught a generation of DSP engineers their craft, either by direct teaching, or by teaching teachers.  They are Alan V. Oppenheim and Richard G. Lyons.  In Oppenheim's case, the Signals and Systems courses that he taught are still available on MIT's web site for free viewing, and his Discrete-Time Signal Processing is still the bible for many engineers on the topic.  Lyons still teaches and writes from Besser Associates, and his "Understanding Digital Signal Processing" book is simple and lucid.

You can't go too far wrong by copying the masters.


Why does a narrower window in a Fourier transform or STFT give poor frequency resolution?

The other answers are good, but I thought that I would try to give a more intuitive/visual answer since I am an intuitive/visual guy.

The picture below is the plot of two tones that are almost the same frequency.  One tone is plotted in red, and the other in blue.

I generated the picture in Matlab with the following code:

tone1 = sin(2*pi*.05 * (0:99));
tone2 = sin(2*pi*.0501 * (0:99));
plot(tone1)
hold on
plot(tone2, 'r')

As you can see from both the picture and the code, the two tones are very close in frequency.  They are starting to separate a bit towards the end of this short window, but they are still very similar, to the point where a DFT could not tell them apart.

If I create the same tones with the same frequencies, only I make the window longer (5000 samples this time instead of 100) we get a much different picture.  Obviously there will be more cycles in each tone...

... but that's not the interesting part.  We see the interesting part when we zoom in at the end of the window.

We see that the two tones are 180 degrees out of phase at the end of the window, which makes them very easy to distinguish.  So how did I know to pick 5000 samples?  The difference in the two tones is $.0002\pi \frac{radians}{sample}$, so $5000$ samples $* .0002\pi \frac{radians}{sample} = \pi$ radians.

As another answer mentioned, the resolution of the DFT is the sample frequency divided by the number of samples (i.e. the window length).  I implicitly made the sample frequency 1, so if we divide 1 by 5000 (the number of samples), we get that the resolution is .0002 Hz.  Our tones actually differ by .0001 Hz, meaning that tone2 will have energy in both tone1's bin and the bin next to it (and a little bit in all the others too, but that's a different story).


Properties of Spectral Transformations - Allocation (decomposition into even and odd part)

Let's say that your signal is composed of two parts: even and odd:

$$s(t)=s_e(t)+s_o(t)$$

We also know following properties of this type of functions:

Even: $f(-x)=f(x)$
Odd:  $f(-x)=-f(x)$

Let's calculate the time inversion of your signal $s(-t)$ and apply above properties:

$$s(-t)=s_e(-t)+s_o(-t)=s_e(t)-s_o(t) $$

So now let's do the trick and add the: $s(t)$ and $s(-t)$:

$$\require{cancel} s(t)+s(-t) = \color{blue}{s_e(t)}+ \cancel{\color{red}{s_o(t)}} + \color{blue}{s_e(t)}-\cancel{\color{red}{s_o(t)}}=\color{blue}{2s_e(t)} $$

Solve it for $\color{blue}{s_e(t)}$, and you will get:

$$\boxed{\color{blue}{s_e(t)}=\dfrac{s(t)+s(-t)}{2}} $$

In the end let's subtract $s(t)$ and $s(-t)$:

$$s(t)-s(-t) = \cancel{\color{blue}{s_e(t)}} + \color{red}{s_o(t)} - \cancel{\color{blue}{s_e(t)}} + \color{red}{s_o(t)}=\color{red}{2s_o(t)} $$

Rearrange and you will get:

$$\boxed{\color{red}{s_o(t)}=\dfrac{s(t)-s(-t)}{2}} $$


rate 1/128 channel code for 5 bits of information

Since you say you want to take the channel conditions into account, one thing that immediately comes to mind is to suggest an adaptive link. Assuming you have access to the channel conditions, you can do two things to adjust the link:

Change the modulation: for example, think about using a higher order scheme like 16 QAM for "clear" channel conditions and using QPSK or BPSK for worse channel conditions.
Change the code rate of the channel code: this can be done by choosing a code with some rate, then depending on the channel you can perform puncturing to remove some parity bits (thus increase rate) when the channel is "clear".


What you see are time-domain samples that correspond to the displacement of a speaker from some rest position as a function of time. For example, the first sample corresponds to the membrane's position at $t=0$, the next represents its position at $t=T$, and so on, assuming that there are $T$ seconds between samples. Common sampling intervals for audio include $T = \frac{1}{8000}$ sec and $T = \frac{1}{44100}$ sec.

In order to turn a stream of samples into real audio, a digital-to-analog converter (DAC) is used first to generate a time-varying analog voltage that corresponds to the samples. This analog voltage would be applied to the input of a power amplifier that drives the speaker. The speaker is designed in such a way that the position of its membrane varies proportionally to the speaker's input voltage. At the end of this chain, then, you have a speaker membrane that moves in accordance with the original audio samples, producing the desired audio via its vibrations.

As you noted, the sample values are normalized to a range of $\pm1$. The resulting loudness of the audio can be affected by many non-digital factors, such as the gain of the amplifier driving the speaker or the DAC's output voltage range. The peak sample values $\pm1$ just correspond to "maximum membrane displacement that the hardware can generate." Therefore, there is no way to encode absolute loudness of the audio in such a stream of samples, hence the use of such a normalization scheme.


If you sample a random process for a specific t, you will get one realization of a random variable. For another t, you get another realization of that random variable. This random variable has its statistics which is almost impossible to learn in real world because not all sample paths are observable. See the brown rectangle in the figure below.

That becomes possible in ergodic processes. An ergodic process is the one where the time average of a sample path is the same as the statistical mean. See the purple line in the figure below, which was the observed sample path such as a particular realization of noise in a communications receiver.


Calculating FIR filter coefficients with the windowing method and Convolution algorithm

draw your desired frequency response.  uniformly sample that frequency response $\frac{N}{2}$ times from DC up to just below Nyquist.  reflect those $\frac{N}{2}-1$ samples (not the DC values) as complex conjugate in the upper half of the frequency response above Nyquist.  do the inverse FFT.  if it's rectangular window, you're done.


Try a samplerate of 48000 instead of 44100.  Or vice versa.


You are correct. If the region of convergence of a right-sided signal (like the one you have) does not include infinity, then the signal is not causal.


How does one calculate a pole-zero plot?

Suppose you are given a system with transfer function

$$H(z)=\frac{(1-3z^{-1})(1-7z^{-1})}{(1-4z^{-1})(1-6z^{-1})} $$

Poles

Poles are the values of $z$ for which the entire function will be infinity or undefined. So, they will be the roots of the denominators, right?

Look here, what values of $z$ will turn the transfer function tend to infinity? Obviously it's $z= 4$ and $z=6$, because if you let $z$ equal 4 or 6, the denominator will be zero, which means the transfer function will tend to infinity.

Zeros

Zeros are the values of z for which the transfer function will be zero. As you have guessed correctly, zeros come from numerator. In this case, zeros are $z= 3$ and $z=7$, cause if you put $z= 3$ or $z=7$, the numerator will be zero, that means the whole transfer function will be zero.

So here poles are $z=4$ and $z=6$, and zeros are $z=3$ and $z=7$.

Pole-Zero Plot

About finding the Pole zero plot, you draw a complex plane. Then you put the values of poles as 'X' marks and zeros as 'O' marks. Here I took the liberty of drawing the pole zero plot of the system:

So, for low pass filter, you find out the transfer function, then the poles and zeros.


Significance of steps involved in flipping process of an image?

Looking at the effect you have achieved; namely a flip along the horizontal and vertical axes, all you need is to apply steps 2,3, and 4.

The following MATLAB / OCTAVE code would also generate the same flipped image in 3 steps as would given by your 5 step procedure:

% I is the image data
If = fft2(I);             % S1 - Take DFT2 of the real image
Gf = conj( If ) ;         % S2 - Take conjugate of DFT of image
g  = real( ifft2(Gf));    % S3 - Invert back into spatial domain

figure,subplot(2,1,1)
imshow(I);
subplot(2,1,2)
imshow(g);

The flipping result can be seen from the following symmetry property of DTFT / DFT :

$$f[n,m] \leftrightarrow F(\omega_1,\omega_2) \implies f^*[-n,-m] \leftrightarrow F^*(\omega_1,\omega_2) $$

For real data $f^*[-n,-m] = f[-n,-m]$.


This was my recent experience where I found it was absolutely necessary to use PID control versus the PI control I would typically use: I was developing a TEC (thermo-electric cooler) oven control circuit that had extremely low frequency poles (sub Hz). With a typical PI control the response time would be limited by the lowest frequency pole. Such that in my case the resulting system with the dominant lower frequency pole would take minutes to settle: this can be estimated using the relationship for a first order system of the 10% to 90% settling time in seconds given by $t_r = 0.35/BW$ where $BW$ is the frequency in Hz of the pole. If a single pole is dominant a first order approximation will be reasonably accurate.

With the use of a PID solution, the achievable rise time can be extended out to the next higher pole, and in my case this resulted in a TEC control of under 10 seconds where the plant (TEC) on its own would take minutes.

In practice what you observe if how the "D" enables a very large overshoot of current (the TEC is current controlled), so as long as the current source can supply the needed current, the temperature is very quickly changed much more during transitions than the current that would be commanded under PI control only. This makes sense since the D is reacting to the differential change of the input.

So in summary I am giving one example where it was justified in using a PID control and what the decision factor was: I was faced with a plant that had two lower frequency poles and wanted to have a response time and settling beyond what would be limited by the lowest frequency of the two. By using PID, the response time and settling was then limited by the higher of the two frequencies resulting in a significantly improved response time.

For your last question, no a PID controller does not need a model of the process as you can experimentally determine the gain coefficients in most practical cases (specifically I outline an approach below for systems dominated by lower frequency poles). The process to experimentally set the gain coefficients is as follows:

Set P, I and D all to be zero and then monitor the step response (step the input and measure the output) while slowly increase P only.

Increase P to the threshold of instability-- this does not mean ringing but a sustained oscillation where the amplitude neither further increases or decreases. Once determined, divide P in half and use that as the gain value for P.  There will be expected offset errors which is why we need I.

With P set as determined above, slowly increase I to the threshold of instability. Then divide this value in half. With a typical pole dominant system such as I have experience with, the rise time will be quite fast but there will still be a significant ring down and here we will see the benefit of adding D instead of simply reducing I to eliminate the ringing (which would also correspondingly increase the rise time).

With P and I set as above, slowly bring up D to achieve the desired step response balancing an overshoot with rise time, and final settling time. You will see that the rise time as set above will not be impacted but the ringing will be eliminated. The ideal transition that has the best compromise will typically have a 8% to 15% overshoot. At this point all values can be adjusted from these starting points to optimize given the actual response of the system.

This is a general approach that works with many real world systems dominated by low frequency poles. This will not achieve good results with a more complex system with many poles and zeros in close proximity or an unstable plant with poles in the right half plane, but it is quick to test and determine the practicality of proceeding with this experimental approach.

However, it is certainly very practical and useful to have a model of the process and be able to simulate the operation. Nyquist's approach of modeling control loops in the frequency domain is very practical for black box systems in that we needn't have a mathematical model of the plant in terms of polynomial transfer functions or pole zero locations to mathematically determine stability using the Nyquist criteria. What we do need is the magnitude and phase of the plant which we can do experimentally (when the plant itself is stable, meaning all poles in the left half plane) by sweeping the plant with a sine wave and measuring the magnitude and phase of the response compared to the input. This is the transfer function. Having that alone is sufficient to design the control loop and determine all stability factors (again when the plant is stable on its own, otherwise things get a lot more complicated which I didn't cover here). Further what is great about this approach is the transfer function needn't be limited to polynomials with pole and zero locations which means exponential factors such as true time delays can be easily incorporated and stabilized. The system may be too complex or cumbersome to try and obtain an accurate mathematical model and here we see that simply having the magnitude and phase of the black box system will suffice. This was the great achievement in Nyquist's approach as to its practical application.

For more on the Nyquist approach to determining stability see Nyquist plot interpretation when curve hits the origin


You're doing the same length as in the transmitter.

You're ignoring the purpose of the cyclic prefix: it's both to allow the "echoes" from the previous symbol to die off (i.e. it contains H). That's why your CP is always a bit larger than the length of channel you're expecting.


Processing a heterogeneous image-collection for CBIR

From what you said in the comments, it seems that CBIR is not really what you need here.
  I am guessing that you would help the historians more by clustering the images, rather than by doing
retrieval.

SIFT and SURF are a good place to start.  You might also try the FREAK descriptor, which you can compute at the
points detected by SIFT or SURF.

Then the important question is what you do with all those features.  For a problem like this you should consider the bag-of-features approach or the pyramid match kernel.

Edit: you may be right to look for a different interest point detector. SIFT and SURF detectors give you centers of blobs (roughly circular patches of uniform intensity).  Since you are looking at man-made gravestones with Hebrew text on them, corners might work better.  A classic approach is the Harris-Laplacian corner detector.  Note that the original Harris detector is not multi-scale.  For this application you definitely need multiple scales.  Other corner detectors such as FAST and AGAST have been proposed recently.  I have heard good things about AGAST, but I have not tried it.


Power Spectrum of the Current of Compact Fluorescent Lamp (CFL) Compared to the AC Input

The simple explanation is that the Compact Fluorescent Lamp (CFL) acts as Non Linear System and hence generates signals with different frequency than the input.

Another thought would be that the sampling channel as a whole has different response to those frequencies and it damps 50 [Hz] / amplify 150 [Hz] (Or both).

You take into consideration the whole chain until you get the digital signal.


Establishing the function of a decaying signal

There is Eureqa : Eureqa (pronounced "eureka") is a software tool for detecting equations and hidden mathematical relationships in your data

EDIT : Now below is what I wrote in first place but reading myself again, I haven't really answered the exact question you've asked ... now that I (think) have (above), I just left it because in  my case Eureqa yielded no good results. (pretty much the same issue as you, I was trying to estimate when a turntable would really stop for my DJ-ing software)

So here it is :

What about a http://en.wikipedia.org/wiki/Linear_regression ? You might be interested in the Estimation category, on the right of the page.
There's the http://en.wikipedia.org/wiki/Zero-crossing_rate, it is quite primitive but yielded some nice results for me, you can find out how many times and how far occurences do occur for your value of 100. (You might need a lot more data than a few samples, though ... but you can oversample your signal)
Last thing, you can (and probably should) apply some filtering to the incoming signal, I did use the Moving Average Filter http://www.dspguide.com/ch15/5.htm. Dependings your needs, you can get rid off some of the noise present to even enhance your estimation.

From this great book : The Scientist and Engineer's Guide to
Digital Signal Processing
By Steven W. Smith, Ph.D. It is free ! Chapters of interest for you : 14, 15, 16, 19, 20, 21

Using Wikipedia and this book, I've been able to track a time-code signal for my DJ-ing software.

(if you don't know what a time-code is : http://en.wikipedia.org/wiki/Vinyl_emulation_software)

Also it's a lot of trial and error, magic numbers sometimes, in my case it was 21 samples (for the MA filter) at 44.1Khz, with this I've been able to track changes in the frequency of about 0.25Hz which was way better than when not using the MA filter before the ZCR.


How to evaluate Spectral Density of DC component vs. AC component during computing PSD?

PSD is the spectral power density and also has the DC component, but in practice finite length of window inevitably produce leaking (because of side-lobes), so to prevent masking of small amplitudes near zero frequency we have to remove DC component.

Functions like pwelch don't remove DC component but it's better to remove it manually and after PSD calculations find it manually.


Find output of digital filter given input and transfer function?

Standard difference equation of a biquad: y[n] = b0*x[n]+b1*x[n-1]+b2*x[n-2]-a1*y[n-1]-a2*y[n-2];

Initialize all x[n], y[n]  to 0 for n < -1

b0,b1,b2 are 0.86 -1.72 0.86 (numerator)

a0,a1,a2 are 1 -1.71 0.73 (denominator)


What do the filter coefficients in a digital filter represent?

We can try a very short introduction:

Every filter represents a Linear Time Invariant System (LTI)
Every Linear Time Invariant System can be completely described by it's transfer function or it's impulse response. The two can be converted into each other by the Fourier Transform
Filter coefficients are derived from impulse response or transfer function
The exact nature of the filter coefficients depends on the algorithm (there are quite a few of those)
In the case of the simplest algorithm, the direct convolution FIR (Finite Impulse Response) filter, the filter coefficients are simply the impulse response of the LTI system.
In most other algorithms the relationship is much more complicated and text book study is indeed required.
The whole subject of LTI systems, transfer functions, Fourier Transforms, amplitude responses, phase responses etc. is probably another text book worth of stuff


Suggestions on how to use FIR model in sentence retrieval task

How can I use the above model where the input signal is a sentence.

makes no sense; your whole system is not linear (which is one of the requirements for your whole convolution approach); for example, if you you play back the words "fox" and "jumped" at the same time, you don't hear "lazy".

Not linear means your system can't be represented by an impulse response.

Your signal model is overly simplistic, too. You might map words to numbers as soon as you know the words, but that's only happening after speech recognition (and might be helpful in understanding logically what has been said), but not before.

So, your whole approach shows you haven't dealt with a lot of speech recognition before – which is OK, but it means that you can't just come up with some approach and assume it makes sense from the start.

I'd recommend reading a book on speech recognition. There's certainly a big wealth of literature out there.


There are some excellent papers on this topic you might want to read.

Sampling-50 years after Shannon, M Unser - Proceedings of the IEEE, 2000

Splines: A perfect fit for signal and image processing, M Unser - IEEE Signal processing magazine, 1999

Interpolation and decimation of digital signals—A tutorial review, RE Crochiere, LR Rabiner - Proceedings of the IEEE, 1981.


Purpose of Extended Kalman Filter with Constant Minimum MSE Matrices

Are there good reasons not to compute these matrices?

Yes.  The computation is expensive, particularly for large systems.  In addition, proving stability, particularly for an extended Kalman filter, can range from hard to impossible.
A KF where you don't update those matrices is a known thing; it's called a "steady-state Kalman filter" and it's based on the notion that for some systems $M_{n|n-1}$ and $M_{n|n}$ will asymptotically approach constant values.  To construct one, you compute* the values of $M_{n|n-1}$ and $M_{n|n}$ in the limit as $n \to \infty$, and use those rather than the continually-updated ones.
Depending on the system, a steady-state Kalman filter can be stable, well-behaved, and much less computationally expensive than a full-on Kalman.  The downside is that it takes longer for the errors to settle out, and in the case of an extended Kalman filter that asymptotic performance may not be as good as it could be because the filter doesn't match the system as well -- but for a lot of systems the provable stability and lower computational burden may be Very Good Things.
* You compute these either by just running the computation forward offline and saving the values, or by using the Ricatti equation to compute them directly.


i will reframe the question a little:

we have

$$\begin{align}
 x(t) &= \sum\limits_{n=1}^N r_n \cos(\omega_n t + \phi_n) \ + \ \epsilon(t) \\
\\
  &= \sum\limits_{n=1}^N a_n \cos(\omega_n t) + b_n \sin(\omega_n t) \ + \ \epsilon(t) \\
\end{align}$$

where

$$\begin{align}
 a_n &= r_n \cos(\phi_n) \\
 b_n &= -r_n \sin(\phi_n) \\
\\
 r_n &= \sqrt{a_n^2 + b_n^2} \\
\\
 \phi_n &= \arg\{a_n - j b_n \} \\
        &= -\operatorname{atan2}(b_n,a_n) \\
\end{align}$$

and $\epsilon(t)$ is the unknown error signal.

first you have to assume that $\epsilon(t)$ contains no components having frequencies of any of the $\omega_n$ in it, otherwise that component will just team up with the component having amplitude $r_n$.

if you know all of the $\omega_n$ in advance, you can calculate all of the $a_n, b_n$ by using an inner product:

$$\begin{align}
a_n &= \frac{\omega_n}{M\pi} \int\limits_{t_0}^{t_0+\tfrac{2M\pi}{\omega_n}} x(t) \cos(\omega_n t) \ dt \\
\\
b_n &= \frac{\omega_n}{M\pi} \int\limits_{t_0}^{t_0+\tfrac{2M\pi}{\omega_n}} x(t) \sin(\omega_n t) \ dt \\
\end{align}$$

$t_0$ can be the most convenient time for you and $M$ should be a large positive integer (and need not be the same $M$ for each frequency component $\omega_n$), but the wider your integration is, the better error rejection you will have.

from the $a_n$ and $b_n$ coefficients, you can calculate $r_n$ and $\phi_n$.


Why is a precoder necessary for DQPSK and what does it accomplish?

Regard the dibits as Gray code representations of the
integers $0,1,2,3$, more specifically,

$$[0, 0] \leftrightarrow 0, ~~ [0, 1] \leftrightarrow 1, ~~ [1, 1] \leftrightarrow 2, ~~ [1, 0] \leftrightarrow 3.$$

Then, the precoding scheme is nothing but differential encoding for QPSK
with

to send = prev - to_encode modulo $4$.

For example the line [ 1 0 ] [ 0 1 ] [ 1 1 ] in the question says that to encode $3 = $ [1 0]  when the previous symbol was $1 = $ [0 1], we send $1-3 = -2 \equiv 2 \bmod 4$ where $2 =$ [1 1].
As stated in the answers and comments
on the previous question,

  Differential encoding for QPSK is not the same as differential
  (binary) encoding of the I and Q bit streams separately.

So the precoding scheme is just differential encoding for QPSK.
As to why it is included in the system, it must be that the
channel characteristics and system requirements are such that
coherent detection of the signal is not thought to be feasible,
and so the receiver is being implemented
is a differentially coherent detector for DQPSK (differentially encoded
QPSK).  Note that coherent detection requires that the receiver
phase lock loop acquire (and track) the incoming carrier in
both frequency and phase (including resolution of
the four-fold phase ambiguity) while for differentially coherent
detection, frequency acquisition and tracking suffices, so that
the receiver is simpler and more robust, though of course
for a given SNR, the BER of the differentially coherent
receiver is higher than the BER of the coherent receiver.


Just get a multi-channel recording interface and some regular analog microphones. For example http://www.sweetwater.com/store/detail/Scarlett18i20/ has plenty of channels and microphone pre-amps.


Taylor series expansion in mean shift tracking

For $y\approx y_0$ you have
$$\begin{align}\sqrt{p_m(y)q_m}&\approx\sqrt{p_m(y_0)q_m}+\frac{p_m(y)-p_m(y_0)}{2\sqrt{p_m(y_0)}}\sqrt{q_m}\\&=\frac12\sqrt{p_m(y_0)q_m}+\frac12 p_m(y)\sqrt{\frac{q_m}{p_m(y_0)}}\tag{1}\end{align}$$
So the function you're approximating is $f(x)=\sqrt{x}$, with $x=p_m(y)$, and you do this in the vicinity of the value $a=p_m(y_0)$.


How to obtain digital IIR filter coefficients for known z-plane zero-pole plot?

This might look trivial, but with the pole-zero diagram alone one can get the transfer function with a scale (gain) ambiguity.

That is, all transfer functions $(Kb,a)$, where $K$ is a constant, have the same pole-zero representation. For this reason, we usually use pole-zero-gain terminology.


The formula you wrote corresponds to the error-free capacity. Namely, Shannon is telling us that $C$ is the maximum (theoretical) limit one can achieve.

Shannon managed to prove that if the information rate on the channel being used is less than $C$, then it is theoretically possible to find a coding system for the signal such that the transmission will have an arbitrarilly small probability of error. Unfortunately, the theorem does not tell us how to find that code, just that it exists.

That's why there is no reference to error rate in the formula for capacity. The formula itself is telling you the maximum information rate you can have with an arbitrary small amount of errors, and that value corresponds to $C$.

Sidenote: remember that the formula is only valid for white Gaussian noise.


How to zoom an image using DCT transform

The following IEEE pdf discusses image interpolation (aka zooming) by DCT. Have a look at it
The main idea remains the same: the DCT spectra is zero padded and converted back into time to obtain the interpolated image which yields the zooming effect.


In theory you just do it the same way. The transfer function of a weighted sum of filters with weights $w_i$ is just
$$H_{total}(z)=\sum_iw_iH_i(z)\tag{1}$$
where $H_i(z)$ are the individual transfer functions.
However, in general the order of the total transfer function is the sum of the orders of the individual filters, unless some filters share the same poles. Consequently, the resulting filter order can become too large for a practical implementation. If the orders of the individual filters "are generally around dozens" then you have a problem, even using floating point arithmetic.


Why is a bandpass filtered analog PM signal still PM and not AM

Conventional AM and narrowband PM signals look quite similar:
$$x_{AM}(t) = A \cos(\omega_c t) + m(t) \cos(\omega_c t) $$
$$x_{NBPM}(t) = A \cos(\omega_c t) - A k_p m(t) \sin(\omega_c t) $$
And based on this, their spectrum also look quite similar. However, there's an important distinction between them:
AM modulation changes (modulates) the amplitude of the carrier wave, but keeps carrier phase unaltered. Whereas, PM modulation, NW or WB, keeps the amplitude of the carrier constant, but alters (or modulates) its phase according to the message signal.
Ther reason for this different behaviour is about the way the sidelobes (spurs) are added to the main carrier. In AM, the sidelobes are added in-phase with the carrier (using two cosines for that), while in NBPM they are added in quadrature (sine added to cosine). This quadrature additon, leaves the amplitude (almost) unaltered but changes the phase-angle of the carrier.
So if you treat NBPM as an AM and try to demodulate it that way, you will get constant output (or nothing after the DC blocking), as the amplitude of the NBPM carrier is not changing...


The coefficients $a_i$ of an autoregressive (AR) model of a signal $x(n)$ enter the model definition in the following way:

$$x(n)=\sum_{i=1}^Na_ix(n-i)+e(n)$$

where $e(n)$ is zero-mean white noise and $N$ is the model order. I.e. the signal is modeled as the output of a linear time-invariant filter with transfer function

$$H(z)=\frac{1}{1-\sum_{i=1}^Na_iz^{-1}}$$

and with input signal $e(n)$. There is no simple interpretation of specific values of the coefficients $a_i$. You have to look at the whole polynomial $1-\sum_{i=1}^Na_iz^{-1}$, e.g. by computing the power spectrum of the AR model of $x(n)$:

$$S_x(z)=\sigma_e^2\left|\frac{1}{1-\sum_{i=1}^Na_iz^{-1}}\right|^2,\quad z=e^{j\theta}$$

where $\sigma_e^2$ is the variance of $e(n)$. So the AR coefficients are the parameters of a specific signal model, and they do not have a simple interpretation comparable to mean or variance of a random process. However, in general they represent a much more detailed description of the signal than mean and variance alone.


MATLAB appears to be applying my filter to the spectrum twice, once on each side

Your filter has a real-valued impulse response, which means it has a two-sided spectrum.  Your filter is only applied once, but that application changes both positive and negative frequencies: -0.3739 to -0.18 and from 0.18 to 0.3729 (modulo the transition region of the filter).  If the filter did not do this, you would have a low pass filter and not a bandpass filter.


It sounds like $W_k$ is just a matrix of 1's and 0's.  When $W_k$ is a 1, then the corresponding range value in $D_k$ is used. If it's a 0, then that range value is not used.

Similarly, $T_k$ should be mostly 1, and 0 when the corresponding $D_k$ values are unreliable.

Does that tally with your understanding?  If not, can you elaborate on your question with more details of what you suspect?


Algorithm that Fills the White Holes Within the Character Image

Based on the blog post - The Paint Bucket in Paint.Net 4.0 (Video) I can tell it uses some edge detection to handle similar colors within a piece wise smooth area.
More information is given in the Paint Bucket Tool documentation.
Usually the way it can be implemented is by defining color metric. How far a color is form another color. If it within the tolerance, it is within the mask.

Morphological operation can fill small holes.
I am not sure, but it won't be too hard to implement something similar.
Probably they convert the RGB data to HSL / HSV Color Model. Then just create a mask based on the distance between the reference point to the surrounding colors.  Now the mask creates a logical (Black and White) mask to be filled.
In case of Paint.NET the algorithm should be something like Flood Fill.

So the steps are:

Convert image to the Color Model of choice (HSV / HSL or even LAB).
Generate a mask using color similarity (Distance Function) on the Color Model. The distance function is parameterized by some Threshold parameter.
Use Flood Fill like algorithm to fill the pixels in the mask using the target color. The mask is basically the stopping condition for the Flood Fill operation.

In MATLAB you have imfill() and Flood Fill.
For something more advanced you may use adaptive morphological functions like Active Contours and Morphological Snakes.


Deciding if the received packet is valid or not valid

CRC is a "definite" indicator of the channel impacting the bits and error in received symbol. So deploying a CRC check would definitely help. Length of CRC would depend on the overall system and data packets.


Why must I take abs of filter frequency response?

Frequency response has two parts: amplitude response and the phase response. Both of these are represented as a complex signal when you get the response from freqz. In order to plot the amplitude response you need to use abs. Otherwise I doubt it only shows you the real part which I think is what you see in the first figure. Note that when dealing with the amplitude response, in most cases you will do a further step and convert it to $\textrm{dB}$ scale.


Change duration/tempo of a .Wav file without introducing pops or pitch shift?

I would suggest the following steps:

identify voices and unvoiced sounds: for this you can use zero-crossing rate and intensity as a basic approach or use more complex methods: cepstrum, harmonicity, others.
the unvoiced sounds can be simply cut, usually without concerns
for the voiced sounds apply PSOLA or WSOLA. For PSOLA there are many implementation and it might be easier. So first you must identify the glotal pulses to be able to identify glotal periods. Then you can "cut" full periods (in a overlap and add sense). Always cut periods in the middle of the phonemes and never near the transition zones (borders)


Transform a data set by exploting the vectorfield

The issue with considering anything as one dimensional in this problem is that they cannot be manipulated as required by the solution.
This can be solved by the optimsation of a function that modifies the $x_i$s towards the direction that minimises some expression of "error". The error here seems to be the distance of each point from the nearest intersection with the known lines.
The minimisation of that error alone is not enough to produce the required result here, but we will come back to this later.
The key "tool" here is the manipulation of the $x_i$s via modifications to the line they reside upon. This makes it a bit clearer why $g(x)$ is not a useful way of seeing things here.
For example, how are we going to rotate a $g(x)$? How are we going to "move it" towards the left or towards the right? If the answers to these questions are "we will just add some offset" or "send it through some rotation matrix", then $g(x)$ cannot stay as it is because it would be unaware of the space it is moving in.
Instead, suppose that all $x_i$s rest on a line running through and parallel to some vector $\vec{a}$. This, being a vector allows us to move it around (to translate the vector, simply add the required $x,y$ offset) and rotate it (to rotate the vector, simply scale its components by the appropriate coefficients). Consequently, we now cannot leave the description of the points as $x_i$ but instead something like $p_i$ taken from the set of points $P$ with each point being a pair of $(x,y)$ coordinates.
To complete the loop we need a way by which to change this line the points are resting on by some factor that is the result of the minimisation. We can do this by combining a translation and a rotation matrix and use the minimisation to inform us about the final $\text{offset}_x, \text{offset}_y, \theta$ parameters required to make the vector $\vec{a}$ coincide with the right points on the "fan" that is created by the $f_i$s.
In terms of practically doing this, you can establish a very simple function for manipulating the $p_i$s as $\text{mod_points}(P,\text{offset}_x, \text{offset}_y, \theta)=p_{i,n-1} \cdot M_{translate}(\text{offset}_x, \text{offset}_y,) \cdot M_{rotate}(\theta)$ where $P$ is the current position of the points, $p \in P$ and $n \in \mathbb{N}$ is the iteration step. Other platforms (e.g. MATLAB, Scilab, Octave, Julia and others) will have similar optimisation functionality of course.
At each iteration ($n$) use this to derive the error function you are trying to minimise. The error function takes the known lines and the positions of the transformed points and returns how far away from optimal is the solution. For example, $\text{fit_error}(f_i, \text{modified_P})$.
This $\text{fit_error}()$ will then be used in place of fun in scipy.optimise.minimise() and your initial $p_i$s (or $x_i$s as you call them now) will be the initial positions x0.
After a number of iterations, minimise() will have moved the line that the $P$s rest upon on the right intersection point with the "fan" that is defined by the (now known as) $x_i$s.
But, if you try to do this, you will notice that the minimisation of the error is not a sufficient condition because it allows the following solution:

Is this a valid solution according to the criterion of minimal error? Yes it is, this line produces points that rest on the known lines. Is it an acceptable solution?....No.
What else should we do to produce a valid and acceptable solution? Add a constraint to the calculation function so that the $P$ set gets closer to the intersections with the $f_i$ lines while at the same time the distances between the points themselves are also preserved. You can model this as two expressions for the error and add them with weight coefficients so that you can tune what sort of priority you want to give over either of the constraints.
Hope this helps.
EDIT:
You have two options regarding preserving the distance of the points:

Calculate the distance from the points, but infer where should the $\vec{a}$ go next. Then, since you do have the line that the points are on and their "inter-distances" are fixed, you can use the orientation of $\vec{a}$ to calculate where each point should be, knowing just the first one. (Subsequent points will be positioned at known distance $d$ from the previous known point towards the direction of $\vec{a}$).

Literally include all distances that need to be preserved in the expression for the error. In this case it would be distances $d(p_n,p_{n+1}), d(p_n,p_{n+2}), d(p_{n+1}, p_{n+2}), d(p_n,p_{n+3}), d(p_{n+1}, p_{n+3}), d(p_{n+2}, p_{n+3})$ where $d()$ is a function that returns the Euclidean distance and $n = 0$ is the index in $P$.


How to choose the frequency grid when using DFT with numpy (FFT)

It really depends on your application and what you want to see and the actual FFT call you make.  Typically if your time series is real valued, the DFT will be conjugate symmetric with respect to negative and positive frequencies and the magnitude information is redundant.  Let assume your data consists of N samples and N is even and you used a complex DFT,  also assume you use C style indexes, i.e. X[0] is the first element in your array. X[0] is going to correspond to the "DC" bin or a frequency of zero.
X[N/2]will correspond to the Nyquist sample frequency.  Most People with ignore X[N/2+ 1] to X[N-1] (or often X[N/2] to X[N])  The thing to remember is that each bin has a conjugate valued twin.  X[0] and X[N/2] are pairs of a sort but don't follow this rule.  X[1] is paired with X[N],  X[2] with X[N-1] , X[3] with X[N-2], and so forth.  One is the positive frequency and the other the corresponding negative frequency.

If you recall your sampling theory, the Discrete Frequency Spectra is periodic, so you can think of the values X[N/2] to X[N] as either the negative frequencies of the spectrum of the cycle centered at DC or the negative frequencies of the cycle centered at the twice the sample (+) frequency. I typically will reorder the upper frequencies so that they are displayed below so the DC bin is the center point of my plot and Matlab has a function called FFTSHIFT that does that, and it also works if N is odd.

If your time series is complex, there are different possibilities.  If x is a Hilbert transpose, the negative frequencies are zero, so unless you want to verify that, it doesn't make much sense to display zeros.

If x is a signal such as an I and Q pair out of a demodulator, using a FFTSHIFT ordering is generally preferable because the frequency shifted carrier will be in the center.

The fact that the upper half of the output are actually the negative frequencies was something that confused me many moons ago.


Support Vector Machine: A non-probabilistic binary linear classifier

It can be said that, it is not based on a generative model which models the generation process of the data and itself a probabilistic model. SVM is based on a minimization problem, which seeks for a hyperplane that optimally separates data points of two clusters given the labels. There is nothing probabilistic in this formulation. You just formulate an optimization problem and solve it via Lagrange multipliers.

I heard (but not studied) that, there are probabilistic SVMs, maybe if you study them, you can see the difference.

Note that, at the fundamentals, authors talk about probabilistic notions, such as risk minimization, however at the macro scale, there is no probabilistic notion.


Defining the SNR or PSNR for color images (3 channel RGB files)

Your PSNR formula has in the denominator a sum of squared distances in a color intensity 1-D space between original and reconstructed pixel values, the sum is taken over all the image pixels. Therefore, you question can be re-formulated as "how one defines a distance in a color 3-D space".
You did give an implicit answer, one of many possible and still valid, in your question: your elaboration on the definition of color space ("each pixel is associated with a 3-tuple (R,G,B)") suggests that you consider the RGB color space. The Euclidian metric can be selected and you have the formula for a color distance between two color vectors F1 and F2
$$
Δ_{RGB}(F_1, F_2) = \sqrt{(R_1-R_2)^2+(G_1-G_2)^2+(B_1-B_2)^2}
$$
and it is true, but not the whole truth.
The practical usefulness of the color distance definition is task specific. From one point of view, it depends on what is the purpose you assign to your image analysis instrumentation. Considering the image impairment by a reconstruction noise, you may want to select a color space adapted to human vision, a perception based color space, as a hue-saturation-intensity (HSI) color space. The metric in HSI used for color distance calculations is
$$
Δ_{RGB}(F_1, F_2) = \sqrt{Δ_I^2+Δ_C^2} \\
Δ_I = |I_1-I_2|; Δ_C=\sqrt{S_1^2+S_2^2-2S_1S_2cosθ} \\
θ = (|H_1-H_2|<π)?(|H_1-H_2|):(2π-|H_1-H_2|)
$$
In the other project, you may want to construct an image from measurement data for presentation in pseudo colors. It will require to work with the proprietary color space that reflects transfer functions of your detectors.
Summing up, your project may require a standard color space, and you simply select the best color space for your task, or it is a sophisticated research work, and to define a PSNR formula for your image processing is a work of its own.
About "color spaces": a search for "distance in color image space" gives as a first entry an article that I readily recommend to you: https://profs.info.uaic.ro/~ancai/DIP/laborator/lab7/color%20spaces%20&%20distances.pdf
In the article, a section 3.4.3 Opponent Color Spaces (page 62, or 26/34 in the pdf) may be what is of special interest for you.


deciphering the frequency response from the input data

Since the signal is most certainly real-valued, you can ignore the FFT values after $f_s/2$ (or sample 300 in your plot) since they will mirror the first ones. With a sampling frequency of 2500Hz, the highest frequency in your signal can be 1250Hz. If you have higher frequencies and don't filter them prior to sampling, you will get aliasing. Look up the sampling theorem for more information, it's quite important when working with sampled signals.

For your actual question I'm not sure what you want to achieve. With an FIR filter you can attenuate certain frequencies. If your transients are strictly not within your desired frequency range, that might work. If they are in the same frequencies as your actual signal, an FIR filter won't help. Is that supposed to be a transient in the first plot where the signal goes from being clean to jiggly?

If so, you could probably just use a lowpass filter which cuts off at 50Hz or similar.


Least-squares digital IIR filter design (with arbitrary responses)

Well, let's decode it:

om=pi*[linspace(0,0.2,20),linspace(0.4,1,60)];

om is a variable that ranges from 0 to pi with two different spacings (one from 0 to 0.2*pi and one from 0.4*pi to pi.

Then, to find D (the desired response), we do:

-1i*om(1:20)*tau

This takes the first twenty entries in om and multiplies them by tau and -i1 = $-\sqrt{-1}$.

That makes the frequency range (D in the above) from 0 to 0.2*pi*tau (i.e. pi since tau is 5).

Checking page 171, we see:

om frequency grid 0<=om<= pi
D complex desired frequency response on the grid om

which means that om has to be between 0 and pi, and D is the desired response.

If you've changed tau, causing D to have a greater phase (more delay).

What I suspect is happening is that the phase (and, therefore, the delay) you've changed it to cannot easily be approximated by a $4^{\rm th}$ order filter (N or M in your script).


How can you perform echo cancellation on features of audio signals?

well, while applying a room impulse response to audio (which is what echoes are) is a linear operation on the time-domain signal, it's generally not on features.
So, you're looking to find a non-linear mapping of a somewhat high-dimensional feature vector space to a "cleaned up" feature vector space with some well-known metric on how well that mapping performs as cleanup....
That sounds exactly like the kind of job that neural networks are good at. And that's not a coincidence: When you look under the mathematical hood of all your ML algorithms, there's the Universal Approximation Theorem that says that a NN that's sufficiently large can represent any non-linear mapping.
So, really, that's a hard problem to do analytically, and it's thoroughly reasonable to assume that "echo removal on MFCC features" can most sensibly be done with a Neural network, most likely even the same that you directly train to classify your audio signals.


How to select the order of LPC to model formants?

Usually you use 2 complex conjugate poles for each spectral peak. By other words, the order of the LPC model should be the double of the spectral peaks you want.
If you have 3 formants then you should use, at least, a 6th order LPC model.
You can try this yourself: record some vowels and then try (for each vowel) to use matlab or python to plot an LPC based spectrum. Change the order of the model to see if it can mimic the original spectrum.


Your question is ambiguous.  Chord, pitch and frequency detection/estimation are very different things.  Bare FFT magnitudes are really bad at guitar pitch modeling.  YIN is a better monophonic pitch estimator.  Polyphonic chord estimation is vastly more difficult, and may or may not even be possible, depending on your input data or time/pitch accuracy requirements (exactly what is your "exact result"?) plus any computational limitations ("real time" latency requirements, CPU performance/cores/GPUs available, etc.)


Speech recognition algorithm in GNU Octave to convert speech to text (and vice versa)

Is there any possibility to implement a speech recognition integrated with google voice that converts speechs in text in GNU Octave?

Yes.
The Google Cloud command line utility returns results in JSON. Octave is working on its JSON support, but until then you can use something like this.
What we are missing now is triggering the whole process from within Octave.
This can be achieved with the Subprocess control functions.
Usually, you can redirect the stdin, stdout streams towards buffers in your application, which would get you the JSON in a local variable within Octave from which then you can call its parser to make it computable. If you cannot, you would have to create an intermediate file and read any data through that. Which might require some timing control too.
This is the general way of going about it even across languages (e.g. Python) and exchanging JSON is probably the better option than just returning some data over a stream and relying on putting it back together on the other side.
But this is also slow which your users (subjects) might find annoying. There is the overhead of starting up the process and in this case, there will be data exchanged over the network too. Besides, the free version of this google service allows a limited number of accesses per month. So, overall, it might be better to check if you can handle that recognition bit locally. You would still probably have to call an external to Octave package (...indeed, most probably Python with the same way outlined here), but the round-trip will be quicker (and cheaper).
Hope this helps.


Independence of Functions of random Variable

It happens that if $X$ and $Y$ are independent then so will their functions $g(X)$ and $h(Y)$ be; but not $g(X,Y)$ and $h(X,Y)$.


with SciPy, how to generate waveform with common pulse shapes?

You have to write your own method to generate the root raised cosine pulse shape. This link might help you: http://en.wikipedia.org/wiki/Raised-cosine_filter

you have to select the parameters for creating one pulse shape and then you can create pulse stream for your base band signal.


What is effect of fractional sample delay effect on FIR Filters?

To answer your primary question, what is the effect of fractional delay- in this case the fractional delay when it would exist (for Type 2 and Type 4 filters only) would be M.5 samples where M is an integer. (So would always be on a half sample increment).

In the picture below the blue line represents the sampled waveform that would pass through the filter, with the blue dots representing the input of the filter, while the output would have samples on the red dots (with the additional M sample delay). This is the effect of a fractional delay.

The presentations are simply showing you the delay of a linear phase filter; in all cases the delay is $N/2$ where N is the order of the filter, which is 1 less than the length L. If the filter has an odd length it will have an even order as in the case of Type 1 and Type 3 filters and the delay will be an integer number of samples. If the filter has an even length it will have an odd order as in the case of Type 2 and Type 3, and therefore the delay will be a half sample more than an integer number of samples.

Consider how this works which also illuminates why these four filter types are linear phase.

The frequency response for a FIR filter is determined by replacing evaluating $H(z)$ as z is swept over the unit circle. ($z=1$ corresponds to "DC" while $z=1$ corresponds to $\pm F_s/2$ where $F_s$ is the sampling rate.) To keep z on the unit circle, $z=e^{j\omega}$ as $\omega$ is swept from $-\pi$ to $\pi$ corresponding to $=F_s/2$ to $F_s/2$.

Therefore, we are evaluating $H(z) = H(e^{j\omega})$

Also note Euler's identities:

$e^{j\phi} + e^{-j\phi} = 2\cos(\phi)$

$e^{j\phi} - e^{-j\phi} = 2j\sin(\phi)$

Note that a real waveform such as $\cos(\omega n)$ has a phase that can only be $0$ or $\pi$ radians (mod $2\pi$).

Similarly the imaginary waveform such as $2j\sin(\omega n)$ has a phase that can only be $\pm pi/2$ (mod $2\pi$).

It is only complex waveforms such as $e^{j\omega n}$ with both real and imaginary components that can have arbitrary phases. And note that the magnitude of $e^{j\omega n}$ is always one, so this expression describes phase only.

From this, the approach is to factor the linear phase frequency response into the form $A(\omega)e^{j\phi(\omega)}$, where $A(\omega)$ is real and only describes the magnitude while $e^{j\phi(\omega)}$ describes the phase:

Given the general form of the z-transform for an FIR filter as

$$H(z) = \Sigma_{n=0}^{N-1}h[n]z^{-n}$$

And the frequency response as evaluated on the unit circle as:

$$H(e^{j\omega}) = \Sigma_{n=0}^{N-1}h[n]e^{-j\omega n}$$

If we factor out of this $e^{-j\omega N/2}$ the result is symmetric terms of $e^{-j\omega n}$, resulting in sine and cosine terms.

Here is a simple example of a 5 tap unity gain filter (coefficients = [1 1 1 1 1]). The is a length = 5 filter that is 4th order (even).

The frequency response of this is given as

$$H(e^{j\omega}) = 1 + e^{-j\omega}+e^{-j2\omega}+e^{-j3\omega}++e^{-j4\omega}$$

as $\omega$ is swept from $0$ to $2\pi$ or alternatively $-\pi$ to $\pi$ depending on how you want to present the frequency axis (corresponding to DC to $F_s$ or $-F_s/2$ to +$F_s/2$.)

If we factor out the linear phase component, we easily get the desired $A(\omega)e^{j\phi(\omega)}$ form since we can replace the remaining symmetric exponential terms with sines or cosines using Euler's identities. In this case, in order for the exponential exponents to be symmetric resulting in real sinusoidal terms, we need to factor out what will result in an integer delay of samples:

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

Notice the symmetric exponential components such that we can now convert them easily using Euler's identities to cosine terms in this case:

$$H(e^{j\omega}) = e^{-2\omega}[(e^{-2\omega}+e^{j2\omega}) +  (e^{-j\omega}+e^{j\omega})]$$

$$= e^{-2\omega}(2\cos(2\omega)+ 2\cos(\omega))$$

So here we see that the term $(2\cos(2\omega)+ 2\cos(\omega))$ contributes nothing to the phase except for $0$ and $\pi$ radian phase changes, so only changes the magnitude in the frequency response. Similarly  the term $e^{-2\omega}$ only contributes to phase as the magnitude is always 1. Notice as the frequency term $\omega$ is changed, the phase is changing in direct proportion, in this case at a negative rate of twice the frequency-- this is consistent with the delay of two samples, and is "linear phase", as the phase is changing linearly with frequency.

We can only extract a single linear phase term and get real sinusoids when the filter coefficients are symmetric or anti-symmetric as in the four linear phase types. Further we see how when we have an odd number of taps with asymmetric coefficients how the center tap must be zero as well as other interesting properties that are revealed. Such as Type 2 filters must always have a zero at z=-1 (low pass) and similarly Type 4 filters must always have a zero at z = 1 (high pass).

It may be helpful to know that time delay is the negative slope (derivative) of phase with respect to frequency--- so if the phase versus frequency response is decreasing linearly versus frequency- this would correspond to all frequency components experiencing the same delay and thus there would be no distortion introduced.

So the reason for showing the fractional delay in those slides (and specifically it will be a 1/2 sample fraction) for the type 2 and type 4 filters is to simply show that this will always be the case for those filter types, and why. So the waveform at the output of a filter that has a half-sample fractional delay will have sample locations that are at the midpoints of the input samples for the underlying waveform that could also be represented in continuous time.


I remember reading a similar question (either here or on stack overflow) abouylt melody extraction. In any event someone did create such a project as a PhD thesis (luckily I bookmarked the site). He even offers his work (for research/personal/non-commercial use) for windows mac and linux http://www.justinsalamon.com/melody-extraction.html#demo if you could somehow integrate this to your application you can use the extracted melody to display the visuals, but still just play the original song from speakers/headphone


To find the filter Co-efficients of the following equation

You're overcomplicating things here. There's no need for sines and cosines and squares. Note that $\omega=0$ corresponds to $z=1$, and $\omega=\pi$ corresponds to $z=-1$. From the definition of the $\mathcal{Z}$-transform you should be able to figure out that the DC term of the transfer function $H(z)$ is given by

$$H(1)=h_1+h_2+h_3\tag{1}$$

and the value at Nyquist ($\omega=\pi$) is

$$H(-1)=h_1-h_2+h_3\tag{2}$$

Equating $(1)$ and $(2)$ with the desired values at these two frequencies gives you two equations for three unknowns. The third equation comes from the linear phase requirement. I'm sure you can take it from here.


SVD is not applicable to time series analysis?

SVD of multi-variate time-series is of course very common.  The point I think the author is making is that SVD cannot meaningfully be directly applied to a univariate time-series, instead it is commonly used indirectly to explore the main eigenvalues of a transition matrix, estimated from observed trajectories of a certain length. see http://en.wikipedia.org/wiki/Singular_spectrum_analysis .  SSA assumes the time-series is stationary and can be modeled as Markov.


integral with vectoric varying coeficient

If you have a closed-form expression for the weighting function, you can use it  directly. Here, the weight is linear, so you can replace coef_vec by xand be ok.

fun = @(x)x.*exp(-x.^2).*log(x).^2;
q = integral(fun,1,10);

If not, an option is to define coef_vec on a larger interval, interpolate it with a suitable closed-form interpolating function, and use this under the integral.

Or, simply, evaluate the integral with a function like trapz. Here is a comparison, with close results

coef_vec=linspace(1,10,100);
fun0 = @(x)exp(-x.^2).*log(x).^2;
fun1 = @(x)x.*exp(-x.^2).*log(x).^2;
q1 = integral(fun1,1,10);

q2 = trapz(coef_vec,coef_vec.*fun0(coef_vec));
disp([q1 q2])

>>    0.0245    0.0245


It feels like any wave at a particular frequency will have a constant phase difference relative to all waves with the same frequency, but from the explanation in the paragraph above that would mean that coherence should be 1 across every frequency for every signal.

I think you're falling into the frequency-domain trap.  The frequency domain is nice, but it was invented to use when thinking in the time domain is harder.  If thinking in the time domain is easier -- don't try to wedge things into the frequency domain, unless you have to for some other part of the problem.

So if you have two signals that are predominantly sinusoidal, and if they're both frequency modulated, then you'll have something like $$s_n(t) = \cos\left( \int \omega_n(t) dt\right)$$ where $\omega_n(t)$ is a random process.  If $\omega_1(t)$ exactly equals $\omega_2(t)$ then your signals will be coherent -- even if $\omega_1$ and $\omega_2$ are time-varying.  If $\omega_1(t)$ and $\omega_2(t)$ have equal long-term averages, but small differences at any given time, then your signals will show some coherence (depending on how you assign numbers to "small differences" and "some coherence").  If $\omega_1(t)$ and $\omega_2(t)$ are unrelated, then your signals will have no coherence.

But note: this is really flipping hard to analyze in the frequency domain and not too bad to think about in the time domain.  So thank M. Fourier for all his hard work, do this analysis in the time domain, and then use the frequency domain the next time you really need it.


What are the ranges of values for each audio feature/descriptor?

Unfortunately, the original standard is not available for free. Furthermore, freely available information on the MPEG-7 descriptors seems to be scarce. Therefore, my recommendation is to get [1] which contains in-depth information about all descriptors, including their computation. Once you understand what each descriptor is supposed to describe, and how it is computed, you'll get an idea about its range of values.

The following resources may also be of interest:

An online computation of a number of MPEG-7 audio descriptors for a wave file

References

[1] H.-G. Kim, N. Moreau, and T. Sikora, MPEG-7 audio and beyond: audio content indexing and retrieval. Chichester, West Sussex, England; Hoboken, NJ, USA: J. Wiley, 2005.


ensure stability against coefficient quantization

your filter is in the form of a standard biquad with 5 independent coefficients ($b_2=0$ in your case):

$$ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{ 1 + a_1 z^{-1} + a_2 z^{-2}} $$

your two poles are the value of $z$ that make the denominator equal to zero

$$ 1 + a_1 z_p^{-1} + a_2 z_p^{-2} = 0 $$

or (multiplying both sides by $z_p^2$),

$$ z_p^2 + a_1 z_p + a_2 = 0 \tag{1} $$

and stability depends entirely on whether your poles (after coefficient quantization) are contained in the unit circle:

$$ |z_p| < 1 $$

or

$$ |z_p|^2 < 1 $$

the solution to Eq. 1 is (quadratic formula):

$$ z_p = -\frac{a_1}{2} \pm \sqrt{\left(\frac{a_1}{2}\right)^2 - a_2} $$

and, if the poles are complex conjugate, then $\left(\frac{a_1}{2}\right)^2 < a_2$ and your complex conjugate poles are

$$ z_p = -\frac{a_1}{2} \pm j \sqrt{a_2 - \left(\frac{a_1}{2}\right)^2} $$

in the case of two real poles, when $\left(\frac{a_1}{2}\right)^2 \ge a_2$, you must guarantee (to ensure stability):

$$ -1 < -\frac{a_1}{2} \pm \sqrt{\left(\frac{a_1}{2}\right)^2 - a_2} < +1 $$

that means both of the following must be guaranteed (after quantization):

$$ -\frac{a_1}{2} + \sqrt{\left(\frac{a_1}{2}\right)^2 - a_2} < +1 $$

$$ -1 < -\frac{a_1}{2} - \sqrt{\left(\frac{a_1}{2}\right)^2 - a_2} $$

what is necessary (but not sufficient) is that $|a_1| < 2$ and what also must be true is (getting rid of the $\sqrt{\cdot}$)

$$ \left(\frac{a_1}{2}\right)^2 - a_2 < \left(1 + \frac{a_1}{2}\right)^2 = 1 + a_1 + \left(\frac{a_1}{2}\right)^2 $$

and

$$ \left(\frac{a_1}{2}\right)^2 - a_2 < \left(1 - \frac{a_1}{2}\right)^2  = 1 - a_1 + \left(\frac{a_1}{2}\right)^2 $$

or

$$ a_2 > -1 - a_1 $$

$$ a_2 > -1 + a_1 $$

so for real poles ($\left(\frac{a_1}{2}\right)^2 \ge a_2$) you must make sure, after quantization that

$$|a_1| < 2$$
and
$$-1 + |a_1| < a_2 \le \left(\frac{a_1}{2}\right)^2$$

for complex conjugate poles ($\left(\frac{a_1}{2}\right)^2 < a_2$), it's much, much easier.  i'll leave this for an exercise for the reader but you can easily show that (for complex conjugate):

$$ |z_p|^2 = a_2 $$

which means

$$ 0 < |z_p|^2 = a_2 < 1 $$

and you still will see that $|a_1| < 2$ is necessary because $|a_1| < 2 a_2 < 2$.  it turns out that $-\frac{a_1}{2}$ is the point on the $\mathfrak{Re}(z)$ axis that is midway between the two poles, whether they be real or complex conjugate.

now, if your digital filter is slewing from one stable complex pole case to another stable complex pole place, then $a_2$ is slewing from one value less than 1 to another value less than one and all filters in between will have $a_2$ less than one and will be stable, at least in a "static" sense.  rapid coefficient slewing can create instabilities in and of itself, but if the slewing stops (at a stable condition), any instability also stops.

so in both real and complex conjugate cases, to have stability in the above biquad you must have

$$ -2 < a_1 < 2 $$

and

$$ 0 < a_2 < 1 $$

for the complex conjugate case $a_1^2 < 4 a_2$

or

$$-1 + |a_1| < a_2 \le \left(\frac{a_1}{2}\right)^2 = \frac{|a_1|^2}{4}$$

for the real pair case $a_1^2 \ge 4 a_2$.

note that, in the real pole pair case, we know that $$-1 + |a_1| <  \frac{|a_1|^2}{4}$$ because $$-1 + |a_1| - \frac{|a_1|^2}{4} = -\frac{1}{4}(|a_1|-2)^2 < 0$$.  so there is always a gap for $a_2$ to fit into.


How would you expect the effect of a denoising model? Do you want it to exceed your ability to observe the noise pattern and remove them, or to replace you and remove the known pattern which you define as noise?

I think most denoising algorithms target the 2nd option. That's the reason why it makes more sense to quantify and compare different algorithms through simulations. The paper assumes the main source of noise comes from speckle, which can be further assumed as Gaussian noise. So they added normal distribution values onto the 'clean' image, and PSNR, SNR, and MSE all get involved with the variation of the image. Such variation is just due to the adding values. They are actually trying to find out which algorithm can reduce such variation to the closest level of 'clean' image. So the comparisons have to be implemented with the known 'clean' image.

Strictly speaking, the 'clean' image in the paper is not clean at all. It is experimental data from ultrasonic B mode scan, and includes the unknown noise. I remember you asked a question in SO very long ago, claiming that you failed to reproduce some of the comparison results from the paper. That's very possible, since as I mentioned above, you can only evaluate your algorithms via simulations, and the 'clean' image should have been generated as a computational phantom. Otherwise you don't even know which part of the image from experiment, other than your added value,  is treated as noise in your algorithm.

As far as I know, SNR can be used as the quantification of experimental data, but it is not calculated that way in your paper. The aim of your algorithm is to enhance the region of interest (organ, tissues in the ultrasound image), making it clearer to the diagnosis. As a result, you'd better choose such regions by yourself, and select other background regions as the  'noise' part in SNR calculation. If you have more than one regions not viewed as background, contrast to noise ratio (CNR) is also suited in your contrast quantification (this is more used in MRI image analysis). Besides, edge spread function (ESF) along with modulation transfer function are often used in optical imaging.

You may also find some books on image analysis useful in image quality evaluation.


How to determine the FFT size of LTE downlink receiver?

You don't need to sample at cell FFT rate. Just shifted the received signal to the desired subcarrier and take a lowpass filtering, then sample at user-specific bandwidth (see attached image). The shifting operation can be implemented by a multiplication with phase ramp in the time domain.

However samping at cell specific FFT rate is optimal because the lowpass filter needs some kind of roll-off.


Audio wave amplitude decreased after FIR filter

As per Peter K. recommendation, I will reformulate my comment as an answer, instead of starting by asking:

I think the fact that the cutoff frequency of your HPF is 17KHz (your tone frequency) can perfectly explain why your signal is being attenuated.

I would recommend you to try to move down the cutoff frequency of the High Pass Filter in order to avoid attenuating your tone.

BR


Why zig zag manner scan is used in DCT for image compression?

You are right! When it comes to quantization, bit allocation or issues like that scanning in zig zag manner seems to be most efficient.

As you said, (and I emphasize) in Natural Signals, lower frequency coefficients contain most of the energy. In a DCT matrix, as you go along a row the frequency is increasing in the first dimension and as you go along a column the frequency is increasing in the second dimension. A high frequency in any dimension probably contains less energy. So for example the
$(2,2)$ entry in an $8\times 8$ DCT matrix $\bf{usually}$ contains more energy than the entry $(1,8)$ in the same matrix.


is root mean square possible to calculate similarity of FFT result?

The Discrete Fourier Transform is a unitary transform (maybe scaled).  As such, it preserves the Frobenius norm of the signal.  That is essentially the same as the Root Mean Square calculation.

In other words:

 rms( x - y ) == C * norm( fft(x) - fft(y) )

for an appropriate constant C depending on your fft implementation.

In Matlab or octave, C=1/nfft:

 n=8192;
 x=randn(n,2)*[1;j];
 y=randn(n,2)*[1;j];
 sqrt( mean( abs(x-y).^2 ) ) - norm(fft(x) - fft(y)) / n

 ans = -5.3291e-15


Determine if system is linear time variant

To be time invariant you would have to show that:
$$y(n-k) = (n-1)x(n-k-1) + (n+1)x(n-k+1)$$
(You shift the $x$ and you get a shift in $y$).
But you can't get rid of the $(n-k-1)$ and $(n-k+1)$.
In general any time you have the time term $n$ by itself it will be time variant unless it cancels out somehow.


White Gaussian noise in the continuous-time case is not what is called a second-order process (meaning $E[X^2(t)]$ is finite) and so, yes, the variance is infinite. Fortunately, we can never observe a white noise process (whether
Gaussian or not) in nature; it is only observable through some kind of device,
e.g. a (BIBO-stable) linear filter with transfer function $H(f)$ in which case
what you get is a stationary Gaussian process with power spectral density $\frac{N_0}{2}|H(f)|^2$ and finite variance
$$\sigma^2  = \int_{-\infty}^\infty \frac{N_0}{2}|H(f)|^2\,\mathrm df.$$

More than what you probably want to know about white Gaussian noise
can be found in the Appendix of this lecture
note
of mine.


How do I find the Hilbert transform (qudrature) of a discrete signal by convolution with a Hilbert kernel?

The short answer: with a finite convolution you can't calculate the Hilbert transform exactly.

The long answer: You will have to make approximations and decide what types and amounts of error you can live with. This will depend heavily on the application, the signals itself and what you are trying to do with the results. This paper may give good guidance on math and potential trade offs http://www.andrewduncan.ws/air/index.html

Per se, the Hilbert Transform is not a smoothing operation. For example the HT of a cosine is a sine, which isn't any smoother.


Do integration/differentiation processes work as simple filters?

Yes, integration and differentiation can be linear filters.
You can start from laplace properties that say:

$ \int_{0}^{t} {x(t)dt} \longrightarrow \frac{X(s)}{s} \\ \frac{d}{dt}x(t) \longrightarrow sX(s) $

So you can find transfer function of integration and differentiation:

$ H_{INT}(s) = \frac{1}{s} \\ H_{DIFF}(s)=s $

You can convert these transfer function in digital IIR filters for example by bilinear transform or others digitization techniques. However you should notice that $ H_{DIFF}(s) $ is not causal, then you must add a pole to the transfer faction far from useful signal frequency, and it begin:

$ H_{DIFF_{causal}}(s)= \frac{s}{\alpha s + 1}$ where $\alpha$ is the time constant of the derivative filter a small real $>0$.

Using bilinear transform $H_{INT}(z)$ begin a trapezoidal integrator, using Euler transofrm $H_{INT}(z)$ begin a rectangular integrator, you can see this difference and the digitization in the MATLAB PID page.

I wrote a simple PID program in C that computes these operations using Euler transform, a PID in closed loop doesn't need to be accurate so Euler works well.

In literature there are a lot of ways to implement derivative filters (IIR and FIR) that workaround in elegant way anti-causal problem, however in a lot of situations you can simply digitize the analog transfer functions to make IIRs and FIRs.


What causes this type of distortion in a constellation I/Q plot?

I don’t think this is necessarily additionally introduced ISI (beyond the ISI of the pulse shaping filter itself, which is zero when there is no timing offset), but may be the result of timing offset error. All the points within a small offset from the ideal sampling location appear to be selected, and the overall pattern of the full constellation showing all samples does not appear to show any distortion.
My theory is we are seeing the results of the overshoot and undershoot you would see in proximity to the ideal sampling location due to the expected response of the pulse shaping filter. What I am describing is better explained with reference to the eye diagrams on the I and Q channel, in particular the case I give below with the higher alpha value where it is quite clear how 4 clusters could result with a small time offset. Notice if the sample is selected at Time = 0 that all the samples will be at the same four I, Q locations, but if we shift the sample time slightly higher there will be 16 possible I Q sample combinations (4 on I and 4 on Q):

Here is the same thing with alpha = 0.25 and the clusters would migrate to a larger cloud:

The OP's case is using alpha =  0.35, so jury is out if this is indeed due to timing offset alone. I recommend reviewing the data with eye diagrams as I have done here to further confirm root cause between residual ISI at the correct sampling location, or timing offset error, or both.


The main issue with the example you gave is that the filter design function cheby1 is returning all NaNs, which isn't going to be a very good filter. The problem is how you're specifying the passband/stopband edge frequencies. This particular function is meant to emulate MATLAB's cheby1 function; the frequencies that you give it should be normalized, such that a value of 1 corresponds to half of the sample rate.

import scipy.signal as signal
fs = whatever_the_sample_rate_of_the_filter_input_is_going_to_be
signal.cheby1(8,0.05,[0.1/(fs/2),10.0/(fs/2)],btype='band',analog=0,output='ba')

I don't have SciPy handy, but that should at least correctly design the filter you want.


How do I find the phase response of a filter from the filter equation?

You need to determine the frequency response of the filter. The phase of this complex function is the filter's phase response. If you have the filter's difference equation, you need to transform it to the $\mathcal{Z}$-transform domain, where a delay of $k$ (with integer $k$) corresponds to a multiplication by $z^{-k}$:

$$Y(z)=c_1X(z)+c_2Y(z)z^{-1}\tag{1}$$

From $(1)$ you can determine the transfer function

$$H(z)=\frac{Y(z)}{X(z)}=\frac{c_1}{1-c_2z^{-1}}\tag{2}$$

The frequency response is the transfer function evaluated at $z=e^{j\omega}$, where $\omega$ is the normalized frequency in radians:

$$\omega=\frac{2\pi f}{f_s}\tag{3}$$

with $f_s$ the sampling frequency, and $f$ the actual frequency in Hz.

From $(2)$, the filter's frequency response is

$$H(e^{j\omega})=\frac{c_1}{1-c_2e^{-j\omega}}=|H(e^{j\omega})|e^{j\phi(\omega)}\tag{3}$$

where the phase response $\phi(\omega)$ is given by

$$\phi(\omega)=\text{atan2}\left(\frac{-c_2\sin(\omega)}{1-c_2\cos(\omega)}\right)\tag{4}$$

For more complicated filters it's usually easier to determine the frequency response (and the phase response) numerically by using an FFT. The (sampled) frequency response is obtained by element-wise division of the FFT of the numerator coefficients of the transfer function by the FFT of the denominator coefficients. The (sampled) phase is then obtained as the argument of this complex-valued vector. In Matlab/Octave, the frequency response can be obtained directly from the filter coefficients using the command freqz.


Is there any practical application for performing a double Fourier transform? ...or an inverse Fourier transform on a time-domain input?

No, taking the Fourier transform twice is equivalent to time inversion (or inversion of whatever dimension you're in). You just get $x(-t)$ times a constant which depends on the type of scaling you use for the Fourier transform.

The inverse Fourier transform applied to a time domain signal just gives the spectrum with frequency inversion. Have a look at this answer for more details.


Covariance matrix $P$ and $Q$ for complex valued data

In general, the equations from your first code block work equally for real and complex values.

Note that the ()' operation is the Hermitian Conjugate (i.e. transpose + conjugate). If you are in the real domain, it only becomes transpose (because conjugation on reals doesn't change them). Hence, you would not need to change anything in the code.

Regarding your covariance matrices: Note that a covariance matrix is always positive (semi-)definite and Hermitian symmetric, i.e. $P^H=P$. This ensures that each single variance on the diagonal is a real value. Also, how would you imagine a e.g. Gaussian with complex variance?

However, there is one thing which might be important for you here: When considering white Gaussian noise in the real domain, you generate random values according to

noise_real = sqrt(sigma2) * randn(N,M);

where sigma2 is the variance of this random variable. If you are in the complex domain, you have to generate them via

noise_complex = sqrt(sigma2/2) * (randn(N,M) + 1j*randn(N,M))

i.e. the variance becomes halfed for real and imaginary domain. Then, when you calculate var(noise_complex) you will get exactly sigma2.

I see from your post that your P and Q denote uncorrelated complex variables. In this case, if they are supposed to have a variance of sigma2 on each element, this is the corresponding code:

P = sigma2*eye(3);


Binarization, and then thinning/skeletonization

Your implementation of the ZS algorithm is wrong. The correct output should be:

Here you can see the output superimposed on the original binary image:

Removing the small branches would require a little more effort and is not exactly off the shelf, but basically you could tag every segment between junction points and end points with its length and prune the ones which are too short. Then you can also take into account the "main" direction and not prune in that direction.


What Does It Mean by Gradient / Divergence / Curl of a 2D Image $ g \left( x, y \right) $?

Original Image

Image Gradient

The gradient shows how quickly the color (or greyscale) changes from pixel to pixel. In very "busy" images, the gradient will change much. In areas of less color change, the gradient will be close to zero.

In the image above, the petals don't change color very much so the gradient is close to zero (or small). At the edges of the petals where the color changes from purple to green, the gradient is large.

Image Divergent

In the sense of image processing, it is the inverse of gradient, this is a indicate of minor gray scale variation.In this aproach i used the channel Red and Green as vector field

Image Curl

Its a like a circular motion of variation of pixels. The curl operation measures the tendency of rotation around the point itself.


First of all, in order to have a Fourier transform, the original signal has to be multiplied by the unit step function:

$$x(t)=e^{-t}u(t)$$

Giving indeed the transform:

$$X(\omega)=\frac{1}{j\omega+1}$$

The way to plot a complex function on the frequency domain is by finding both its amplitude and its phase, and drawing one graph for each.

For the amplitude case, find $|X(\omega)|$, which is:

$$X(\omega)=\frac{1}{\sqrt{\omega^2+1}}$$


How to add a frequency offset to a buffer of MPSK signal for simulation purposes?

Multiply the time domain samples by a time-varying complex exponential to effect a (cyclic) shift in the frequency domain.
$$y[n] = x[n]\cdot e^{j\pi \frac{f_{offset}}{F_s/2}n}$$
Where $f_{offset} \in \left(-\dfrac{F_s}{2}, \dfrac{F_s}{2}\right)$.
This is a result of the Shift Theorem for the DFT.


What exactly is the effect of a Hann Window on the FFT output?

See the following: http://www.dtic.mil/dtic/tr/fulltext/u2/a381695.pdf  especially p.631 Figs. 17+18.  As noted by the author, a 4 to 1 overlap using a Hann or Hamming window allows for proper reconstruction.

Also see: F. J. Harris, “On the use of Merged, Overlapped and Windowed FFT’s to Generate Synthetic Time Series Data with a Specified Power Spectrum,” Proc. IEEE 16th Asilomar Conf., 1982, pp. 316-321.


What you're asking for is:

  given a periodic signal (e.g. square wave), where are the discrete frequency components (i.e. harmonics)?

You'll find that the harmonics you're looking for are the coefficients of the Fourier series expansion of your signal.

So, if you have a signal that has a known Fourier series expansion, or can be constructed from such signals, it's easy.

But how to find the coefficients in general?

This is answered exactly by the Fourier transform.

So, the relevant formulas are, given a signal $s(t)$ with period $P$

the coefficients $c_n$ of the Fourier series expansion $s(t) = \sum\limits_n c_n e^{j2\pi n\frac tP}$: $$c_n=\frac{1}{P}\int_{t_0}\limits^{t_0+P} s(t)\cdot e^{-j \tfrac{2\pi nt}{P}}\ dt$$ describing the $n$th harmonic. You just plug in your square wave for $s(t)$ and get out the coefficients of your harmonics.

and in the general, you don't get discrete harmonics, but a continuous spectrum, and thus, instead of a discrete Fourier series the Fourier transform.

Not posting info on the Fourier Transform here; it's pretty much common knowledge, and you wouldn't benefit – your profile indicates that you're studying physics, and thus, you've already met the Fourier transform, for example when describing electrons in periodic crystals.


Matlab : Exploring the effects of co-efficient quantization

For each of the vectors "a" and "b", you need to find the element with the largest absolute magnitude. (I'm not sure of the specific Matlab code for this; probably something like max (abs (a)).)

Then, you need to pick the smallest number of the form 2N (N is an integer) that is larger than the value determined in the previous step.

From there, the rest of the instructions should be straightforward.


Conceptual Questions on Colored Noise Process

Once you talk about the Spectrum of noise / process you implicitly says it is stationary in the wide sense.
What does it mean to have a signal with uncorrelated samples? Do you understand it means you can't using linear predictor, to have any information from all past samples ion the current one? Usually this is not how signals behave. This is exactly why correlated noise is hard. It behaves like a signal hence it is harder to separate them apart. One way to deal with the model of Signal + Colored Noise is to apply whitening transform on the data, apply White Noise methods and then apply the inverse transform.
Any process where the probability for a specific value depends on the index value.
Here are some specific examples (From Investopedia - An Introduction to Stationary and Non-Stationary Processes):

Noise is used when dealing with signals. We call noise to the non deterministic added (Usually, we have multiplicative noise as well) term of the signal. Since we have some index (Time / Spatial, etc...) for the signal we have it for the noise as well which means it is a process.


Estimating scaling and translation between two consecutive temporal images

I feel there is no the best method to solve this, but if you do not have some specific feature point to track, I suggest using optical flow a below:

compute optical flow of two frames,
find the dominate global optical vector, may be through computing average or median of optical vectors (you do not want small noisy motions)
Select at least 2 vectors parallel and similar to the global vector and fit a rotation*scale transformation matrix.

In case you have some feature points, instead of computing optical flow for all pixels simply use sparse optical flow OpenCV tutorial on Optical flow

This link might also be helpful.


You have a system with the following input-output relation:

$$y(t)=\int_{-t}^{\infty}x(-3\tau)d\tau\tag{1}$$

In order to check whether the system is time-invariant or not, we need to compare the shifted output with the output resulting from a shifted input. The shifted output is

$$y(t-T)=\int_{-(t-T)}^{\infty}x(-3\tau)d\tau\tag{2}$$

Shifting the input means applying an input signal $x_T(t)=x(t-T)$. Note that $x_T(-3\tau)=x(-3\tau-T)$, so the response to the shifted input is

$$y_T(t)=\int_{-t}^{\infty}x(-3\tau-T)d\tau\tag{3}$$

The system is time-invariant if $(3)$ equals $(2)$. So we try to make the integrand in $(3)$ look like the one in $(2)$ by substituting $-3\tau-T$ by $-3\zeta$, which gives

$$y_T(t)=\int_{-(t-T/3)}^{\infty}x(-3\zeta)d\zeta\tag{4}$$

Now $(2)$ and $(4)$ have the same integrand, but we see that the lower integration limits are different. Consequently, $y_T(t)\neq y(t-T)$, hence the system is not time-invariant.


Image Interpolation Using the Yule Walker Equations

The main problem is Auto Regressive (AR) / IIR models are based on the idea of something before and after.
While for time series it is well defined, how will you address that in image for a pixel in the middle of the image?

I addressed some ideas on that in my answer to the question Denoising an Image Using Kalman Filter.

Another idea is to work across rows and then across columns.
Then you can build an AR model based on the samples with unit "Time Parameter" and then ask what should happen for $ n / 2 $ or so.
You can even take it farther, and in order to use the 2D information take $ k $ rows / columns at a time and build a vectorized AR model.
Then it means that for each pixel you'll get $ k $ estimations which you could average or chose any other idea to fuse.

Actually it could be nice idea to have a try...


Calculating phase with noninteger periods of a sinusoid

Do an fftshift (rotate the input vector by N/2) before doing the FFT.  Or for a phase measurement at a fixed point in time, start sampling the data earlier so the N/2 time sample in the FFT input vendor is at the time point where you want the measurement, and then flip the sign of the phase of every other bin of the FFT result (the odd numbered ones).  This will move the phase reference point to the middle of the data window, where an input sinusoid will not have a discontinuity in phase that changes with frequency.  This also allows interpolation of phase between FFT result bins.

The evenness/oddness ratio, and thus the FFT phase, does not change with frequency when a sinusoid's phase is referenced to (for example both generated at and measured with respect to) the center of the FFT aperture.


How does shift and scaling inside of a function affect its Fourier Transform?

For one-dimensional Fourier transforms,
$$\begin{align*}
x(t) &\leftrightarrow X(f)\\
x(t-a) &\leftrightarrow X(f)\exp(-j2\pi fa)\\
x(bt) &\leftrightarrow \frac{1}{|b|}X\left(\frac{f}{b}\right)
\end{align*}$$
all of which are proved by change of variables in the definition of the
Fourier transform
$$X(f) = \int_{-\infty}^{\infty} x(t)\exp(-j2\pi ft) \mathrm dt.$$
It is good exercise to work out the details of each of these results
at least once in one's life, and once one has done so, the two-dimensional case
is as easy as rolling off a log.  If you always look up the answers in a table
of Fourier transform properties, the right result will never come to mind just when
you need it most.

More generally, it is a good idea to remember these results in graphic
terms.  If $x(t)$ is a signal with Fourier transform $X(f)$, and $b > 1$,
say $b = 2$,
then $y(t) = x(2t)$ is a time-compressed form of $x(t)$: while $y(0) = x(0)$,
$y(3) = x(6)$, that is, $6$ seconds worth of $x(t)$ have been squeezed down
into $3$ seconds worth of $y(t)$.  Thus, the frequencies have doubled, and so
$Y(f)$ should be something like $X(f/2)$ (but see below).  More generally,

  Compressing the time axis by a factor of $b$ expands
  the frequency axis by a factor of $b$.

But that is not the complete answer.  One of the more boring results in Fourier
transform that everyone tends to ignore is
$$y(0) = \int_{-\infty}^{\infty} Y(f) \mathrm df.$$
But,
$$y(0) = x(0) = \int_{-\infty}^{\infty} X(f) \mathrm df$$
and since we have stretched the horizontal axis by a factor of $b$
in going from $X(f)$ to $X(f/b)$, we better reduce the magnitude of each
frequency component by a factor of $b$ so that the two integrals (= area
under the curve) work out to have the same value, that is,
$$Y(f) = \frac{1}{b}X\left(\frac{f}{b}\right), ~~ b > 0.$$
The result for $b < 0$ can be worked out similarly.


Upsampling vs downsampling. Which to use when?

You should use the one you need for your problem, when you know which components of your signal are of interest to you.
Let's say you have in your electronic editing an ADC digitizing 40M samples per second to study a heart rate of 70 beats per minute, you are very likely to work with useless information, that's why it will be better to down-sample your signal.
In another case, it can be very useful to oversample, some algorithms to compute an FFT work well with sample sets of power 2, it may be necessary to oversample your signal to apply a fast FFT.
That is why you should know the characteristics of your signal of interest before applying any of these functions.


Is it possible to specify number of taps and cutoff frequency independently?

If you talk about a moving average filter then you only have a single degree of freedom, namely the filter length (i.e., the number of taps). Consequently, given the number of taps, the magnitude of the frequency response is fixed and the cut-off frequency is given.
Note that you could also define the delay of the moving average filter as another degree of freedom, but that wouldn't change the filter's magnitude response and its cut-off frequency.
For general FIR filters with $N$ taps you have $N$ degrees of freedom, so you can clearly choose the cut-off frequency independently of the number of taps. But this is not the case for a moving average filter.


Denoising a signal using eigendecomposition

So I'll get this answer started by saying again that what you are describing reminds me a lot of Singular Spectrum Analysis.  The process is pretty much exactly how you describe: you have some observed signal, $Y[n]$, which is some sort of mixture, and we take $Y[n]$ and create a Hankel matrix from it using some desired frame length M (this will effectively determine your sub-space resolution).

Now that we have our Hankel matrix, we can either compute $R_Y$ directly and then do an eigendecomposition, or we could go with a Singular Value Decomposition of the trajectory matrix.  You'll find them to essentially be equivalent if we use the left singular vectors; either method is perfectly valid.

So once we get our matrix of eigenvalues, $U$, we need to project them using your typical Karhunen-Loeve method to get what we'll call principal components (this is sort of classic SSA notation, I dislike the naming but I digress):

$P = Y^H U$.

Now that we're here, we can use the singular values to determine amplitude of the signals themselves, and the principal components essentially determine the "parts" of the signal, i.e. $X[n]$ and $\epsilon[n]$.

Let's say we have N principal components, and we know that we only "need" the first one: from our principal components, we perform what is typically called Eigentriple grouping via the following multiplication:

$C = UP^H$

This new matrix, C, is what we call the reconstructive component(s) of $Y$, and is based on which principal components we selected.  We're not quite done yet, becuase $C$ is the same dimension as our Hankel matrix, $Y$, so do get back to our signal, we'll need an additional step.  To reconstruct the series itself, we'll perform what's called diagonal averaging (a bit of a misnomer I think).  The goal here is that we'll finally map back to a single signal, let's call it $Y_out[n]$.  It's in this step that we'll be able to utilize that one-to-one relationship of the Hankel to the series to extract out signals.  This step is a bit long winded, so I've omitted it for brevity, but you can consult Wikipedia or any of the many papers/texts on SSA for more information.

So now that we've briefly gone through it, you might find yourself asking "when will this method actually work".  SSA, like anything, is a tool and it has its limits.  From a signal processing standpoint, it's great at determining individual signals which are fairly narrowband (i.e. they don't cover the entire spectrum) and are "separable".  In the context of SSA, "separable" means that the signals are non-overlapping in frequency.  For a simple example, let's say I have two sinusoids, one at 50Hz, and another at 500Hz.  SSA will be able to easily separate the two signals, and we will be able to isolate the signals pretty easily; we would call these signals separable. However, let's say I change the second signal to be 51Hz instead of 500Hz.  In this case, SSA will struggle, as the signals are very, very close to one another; we would call these signals non-separable.

How do we determine which signals we care about? That's the difficulty with this method.  If you have some appreciable signal-to-noise ratio (SNR), you can perhaps use that to your advantage and identify which signals are of interest and which signals are just purely noise by inspecting the eigenspectrum; true signals as represented by the principal components will likely be associated with large eigenvalues, whereas noise will be relatively low power, and spread.

So, in short, if you have some time series observation, you can use SSA to attempt to discern individual signals/trends.  It's a really powerful tool, and as long as the signals are separable, you should have some reasonable success.  Hope that helps!