Measure the SNR of a signal in the frequency domain?

It can be a little tricky.  If your sine wave happens to fall at an FFT bin center, things are a little easier.  If your sine wave happens to not fall at a bin center, you have to consider spectral leakage.  You'll also need to consider how your window function affects your signal.

A simple technique to estimate the signal power would be to sum the squared magnitudes of the FFT bins with most of the signal's energy (i.e. the largest bins).  To get the noise energy, sum the squared magnitude of all of the other bins.  This will get you close.


Proof of the convolution property of Fourier Series in continuous time

HINT:

Show that with $\Omega_0=2\pi/T$

$$\int_T e^{j\Omega_0(n-p)\tau} d\tau=\begin{cases}T,&\quad n=p\\0,&\quad\text{otherwise}\end{cases}$$

So the sum over $p$ reduces to a single element $Td_ne^{j\Omega_0 nt}$, and the result follows.


What is the Fourier Transform of an integral

When integrating functions with dirac delta you have to remember that
$$\int \delta(x - a) f(x) dx = f(a)$$ so, it works as  change of variables
and for the 2D dirac, integrating over an axis gives the 1D delta
$$ \int \delta(x - a, y) f(x) dx = \delta(y) f(a, x)$$
$$\begin{eqnarray}g(x,y)
  &=& \int_{0}^{T} f(x,y) \delta(x - \nu_x t, y) dt \\
  &=&\delta(y) f(\nu_x t, y)
\end{eqnarray}$$
This does not depend on $x$ or $u$, it is expressed in terms of $f(\nu_x t, 0)$, if you want to express in terms of $F(u, v)$ you must calculate the inverse fourier tranform for these for this point.
$$\begin{eqnarray}G(u, v)
&=& \int \left(\int \delta(y) f(\nu_x t, y) e^{-2i\pi i vy} dy \right) e^{-2i\pi i ux} dx \\
&=& \int \left( f(\nu_x t, 0) \right) e^{-2i\pi i ux} dx \\
&=& \delta(u) f(\nu_x t, 0) \\
&=& \delta(u) \int \left( \int F(u,v) e^{-2i\pi u v_x t} du \right) dv
\end{eqnarray}$$


The system is not time invariant as can be seen by the fundamental test procudere:

Given the system with an input/output relation:
$$y(t)=\mathcal{T} \{ x(t) \} = \int_{-\infty}^{x(t)} x(\tau) d\tau$$

First check the new output for a shifted input $x(t-d)$ as
$$y_d(t) = \mathcal{T} \{ x(t-d) \} = \int_{-\infty}^{x(t-d)} x(\tau-d) d\tau $$
$$y_d(t) = \mathcal{T} \{ x(t-d) \} = \int_{-\infty}^{x(t-d)-d} x(\tau) d\tau $$

And then also check if the output $y(t)$ shifted by $d$ yields th same formula?
We know that output $y(t)$ is given in terms of input $x(t)$ with the formula
$$y(t) = \int_{-\infty}^{x(t)} x(\tau) d\tau$$

Hence the shifted output $y(t-d)$ is obtained by simply replacing every $t$ variable with $t-d$ in that formula which yields:
$$y(t-d) = \int_{-\infty}^{x(t-d)} x(\tau) d\tau $$

And now looking at the formulas for both $y(t-d)$ and $\mathcal{T} \{ x(t-d) \}$ given above, we see that since $x(t-d)-d \neq x(t-d)$ for all $d$ then $y(t-d) \neq y_d(t)$ and we conclude that the system is not time invariant.


Just "eyeballing" your signal, it looks like the interesting things are all down below 1Hz.  The peaks in your signal browser diagramm seem to be spaced at every 5 to 10 Seconds, which would be 0.1 to 0.2 Hz.  Your spectrum plot goes from zero to 25 Hz, and doesn't show anything really interesting.

Try zooming in on the spectrum and see what lies below 1 Hz.  Yo may well find that the signal itself lives down in the very low frequencies, and that separating the noise from the signal is as simple as using a low pass filter.


Least Squares and Auto and Cross Correlation

Cross-correlation and auto-correlation are not directly linked to least squares.

Least squares can create a line of best fit for a set of points.

Cross correlation is a measure of the similarity of one signal to another (as a function of time lag to one signal). The discrete (digital) cross-correlation is an array where each element tells you how much the signals correlate with a certain amount of time lag. This operation is performed naively with a sliding window, or faster using Fourier, or other transforms. xcorr(a,b) = IFFT(FFT(a)* conj(FFT(b)))

Auto correlation is a measure of the similarity of a waveform to itself as a function of time.

There are normalised variants of cross-correlation and auto-correlation. These are created by measuring the standard deviation of each signal and dividing the output by the multiple of both. This enables you to classify the level of correlation much more easily. 0 implies no correlation and 1 implies identical signals.


Is there any simple algorithm to convert a sound source in 3d to volume and panning

inverse-square law for distance.  and Blumlein stereo for pan.


Given the Graph of a Fourier Series $\sum c_k e^{2\pi ikx}$ Find the Graphs of $\sum c_{3k} e^{2\pi ikx}$ and $\sum (c_k)^2 e^{2\pi ikx}$

This doesn't require any complicated computations! Look up some tables of identities about the Fourier Series. Which operation would you have to apply to $f$ so that its Fourier series would become $\sum_{-\infty}^{\infty}(c_k)^2 e^{2\pi ikx}$?

For example, we know that derivating in the time domain is equivalent to multiplying the Fourier series coefficients by $i k$ ; so we know that the graph of the Fourier series $\sum_{-\infty}^{\infty} ik (c_k) e^{2\pi ikx}$ is the graph of the derivative of $f$.


What does the uncertainty principle say about recursive filters?

The uncertainty principle basically says the same thing for IIR filters: "the steeper the filter, the longer the Impulse Response". The only difference between FIR and IIR is that you have to use a somewhat different definition of "filter length".

While all IIR filters are technically infinite, they have very different decay rates. You can define a quantity that basically describes "energy decay over time". For most IIR filters this is an exponential decay (often dominated by the pole-pair closest to the unit circle) and the decay rate is directly related to the filter bandwidth.

An interesting quantitative approach to determining decay rates can be found in the area of reverb time calculation, using "time reversed impulse response integration".


Why doesn't the sampling theorem work in this case?

The formula given is

$$ F_s > 2F_{\mathrm{max}}$$

not

$$ F_s \ge 2F_{\mathrm{max}}$$

otherwise it is possible that you could be sampling exactly on the zero crossing points. In practice you would choose sometime a bit higher. For example, CD audio uses 44100 Hz to represent signals up to 20000 Hz.


Problem using pitch shifting with TD-PSOLA and formant preservation

I believe you're correct, it's an issue with the overlap-add process. The code doesn't correctly calculate the triangular windows for lower pitches.
In short, the windows are too big.
It produces incorrect results for any pitch ratio below 1, such as 0.5 (one octave down).
The diagram below explains the issue as I understand it:

The big windows introduce so many artifacts (additional peaks) that the result sounds almost exactly the same, apart from the "speaking into fan" effect.
To get around this, it would be easiest to add an if statement which would use the existing windows only when the pitch ratio is equal or greater than 1, otherwise it should use the 'ideal windows'.
However, I'm not sure to implement the 'ideal windows' in Python. The existing algorithm uses newPeaks to calculate windows. For lower pitches it must use the old peaks array. However, this requires new math which I'm not sure how to implement exactly.


Dominant Orientation using Fourier Power Spectrum

Subtract from the 8x8 image the mean value of its pixels, which should not affect the dominant orientation. This will remove artifacts due to windowing in the next step.
Window the image with a circularly symmetric window function that eliminates artifacts aligned with the image axes. Make sure you have chosen a window function that is sufficiently smooth to be sampled with a 8x8 grid without frequency domain aliasing.
Pad the image to 32x32 with zeros to increase the image size and thus the discrete Fourier transform (DFT) resolution. This makes the DFT result smoother and thus easier to interpolate.
Calculate for each angle an integral of the squared magnitude of the DFT result over a line corresponding to the angle. All angles must be treated equally; use same length lines. 4x4 piece-wise bicubic polynomial interpolation should suffice if you approximate the integral by a sum of point samples at the axis-aligned resolution of the DFT.

If you want to make the computation faster you can take advantage of that the integral will change quite smoothly as function of the angle, so you can first sample at the coarse resolution of at least $\text{atan}(1/8)\times 360\unicode{xb0}/(2\pi) \approx 7\unicode{xb0}$ to find where the highest peak might be, and then sample more finely around that angle.


You have two methods. The second one is right, because it's actually upsampling: inserting zeros, and then suppressing images (what you call aliases here is actualyl called images) through the desired filter.
The first one is actually the same, but instead of the RRC, you use a rectangular filter: When you think about it, holding the last value is like inserting zeros, and then convolving with a 1-amplitude rectangle that is as long as the number of zeros you've inserted.


Fourier Transform gives unexpected results: signal reversal and conjugation

Try plotting the signal and the straight reversed signal in points rather than lines and you will see that, unlike with reversal via FFT, none of the values stay the same throughout the reversal. The straight reversal flips the 128-long data around index 63.5 or equivalently around index -0.5 if you think of the data as periodic with a period of 128, rather than around index 0 (or 64) as in reversal via FFT.

According to user David Hoffman's tests both

roll(irfft(conj(rfft(y))*rfft(y)), -1)

and

irfft(rfft(y[::-1])*rfft(y))

work for autocorrelation, roll being numpy circular shift.


Converting floating point number to Q-notation fixed point number in C/Java?


What is meant by "information encoded in time domain"?

The simplest (imho) explanation is this.
Consider a time-domain signal which has narrow pulses or sudden jumps in value or on/off switches at certain (unknown!) time instants. Such a signal is said to carry its information in those specific time instants at which such sudden changes occur.
A PWM (pulse width modulation) signal is an excellent example for such class of signals.
If you consider the frequency domain reflection of those sudden changes in the signal amplitudes, you will see that they have wide band spectrums. They are local in time, therefore spread in frequency.
When you filter such a signal with a typical LPF (low pas filter) or alike, which discards parts of the spectrum away and passes only a subset of it, you may effectively be throwing away part of the information that's carried by the signal too.
For example, the specific instants at which the PWM wave switches on/off becomes blurry, and hence the resolution of the information is reduced.
The author seems to mention this issue with filtering of time-domain encoded information...


What do you call the random Gaussian vectors in compressed sensing?

The collection of scalar products is often called measurement. The collection of vectors (each being an atom) can be called measurement dictionaries or sensing dictionaries. This is not specific to random atoms.


As mentioned in the comments, your question is pretty confusing and confused, so I might not be able to answer all your smaller questions, but I'll try to shed some light on the topic. This should help you understand what's going on, and hopefully this understanding will help you answer all further questions by yourself.

In the following I assume two's complement representation. Note that if $M$ is the word length, and $n$ is the number of fractional bits, a number $x$ is represented as

$$x=K\cdot 2^{-n}\tag{1}$$

where $K$ is an integer is the range $[-2^{M-1},2^{M-1}-1]$. Now if you want to change the number of fractional bits to $m$, you have

$$x=L\cdot 2^{-m}\tag{2}$$

Combining (1) and (2) gives

$$L=K\cdot 2^{m-n}\tag{3}$$

If $m>n$, $K$ needs to be left-shifted by $m-n$ bits to obtain $L$. This means that you lose range (i.e. the maximum number that can be represented becomes smaller), so there is a possibility of overflow. If $m<n$, $L$ is obtained from $K$ by right-shifting $n-m$ bits, but for two's complement you need to replicate the sign bit, i.e. if the original sign bit was set, the new right-shifted number must also have the sign bit set. For $m<n$ you can lose precision because the resolution is reduced.

An example makes things easier to see. Let $M=6$, i.e. the integer $K$ is in the range $[-32,31]$. Let the number of fractional bits be $n=3$. E.g.

$$x=\texttt{101.011}=(-2^{5}+2^{3}+2^{1}+2^{0})\cdot 2^{-3}=-21\cdot 2^{-3}=-2.625$$

If we want to change the number of fractional bits to $m=2$, you need to do an arithmetic right-shift, resulting in

$$y=\texttt{1101.01}=-11\cdot 2^{-2}=-2.75\neq x$$

i.e. precision is lost. If instead we want $m=4$ fractional bits we need to do a left-shift resulting in overflow:

$$y=\texttt{01.0110}=22\cdot 2^{-4}=1.375\neq x$$

This is clear because by increasing resolution we have decreased the range of possible numbers to $[-32,31]\cdot 2^{-4}=[-2,1.9375]$, which makes it impossible to represent the original value $x=-2.625$.

There are of course examples where you would neither lose precision nor would you get overflow. E.g. (again with $M=6$ and $n=3$)

$$x=\texttt{001.010}=10\cdot 2^{-3}=1.25$$

For $m=2$ we get

$$y=\texttt{0001.01}=5\cdot 2^{-2}=1.25 = x$$

i.e. no loss in precision, and for $m=4$ we get

$$y=\texttt{01.0100}=20\cdot 2^{-4}=1.25=x$$

i.e. no overflow.

I hope that this addressed at least some of your questions and that things have become a bit clearer. If you haven't done it yet, you should check out the basics here, here, and here.


Improving the performance of the Hough transform

Are you using openCV?
Anyway, if you are running an edge-detection, and then the hough transform on the edges, it might help to smooth/blur the edges with some bluring functions. OpenCV got some different ones for that.

In particular a boxFilter with a 2x2 sized kernel, anchor=Point(-1,-1), and normalize=false.
It will hopefully make all the edges "thicker", filling in gaps and inconsitencies in the process.


MIT 6.003 HW#8 Problem 4 - Fourier Coefficients of Triangle Wave

The symbol $\underset{T}{*}$ in the solution refers to ($T$-) periodic convolution, which is defined as
$$y(t)=\int_{0}^Tx_1(\tau)x_2(t-\tau)d\tau\tag{1}$$
where $x_1(t)$ and $x_2(t)$ are $T$-periodic signals. If $a_k$ and $b_k$ are the Fourier coefficients of $x_1(t)$ and $x_2(t)$, respectively, then the Fourier coefficients $c_k$ of $y(t)$ are given by
$$c_k=Ta_kb_k\tag{2}$$
Your approach to the problem is valid for $k\neq 0$, and your result is actually pretty close because
$$\frac{1-\cos\big(\frac{\pi k}{2}\big)}{4\pi^2 k^2}=\frac{2\sin^2\big(\frac{\pi k}{4}\big)}{4\pi^2 k^2}=\frac{\sin^2\big(\frac{\pi k}{4}\big)}{2\pi^2 k^2}\tag{3}$$
so there's just some scaling problem that I'm sure you'll be able to fix (check your definition of $f(x)$!). Note that you need to compute the coefficient for $k=0$ separately.


Numerical higher order derivatives and time axis

When considering finite differences there are three of them: Forward, backward and central.

Neither of them is "shorter" than the main signal. The zero is not exactly zero, it should be $\emptyset$. That is, "I don't know".

In backward difference, $y[n] = x[n] - x[n-1]$, but at $n=0$, $n-1$ is not even defined and the same applies for $y[0]$. We usually do $y[0]=0$, (which implies that $x[-1]=x[0]$).

Notice here that we wrote $y[n]$. In finite differences, in the discrete world, there is no $x[0.5]$. It is either $0$ or $1$. Central difference might appear to have a $\frac{1}{2}n$ term but that could be seen as a "phase shift" with respect to the sampling of the signal. The relative temporal difference between two successive samples is still one sampling period.

And this is the interesting bit here. The finite difference is the result of something that happens to two samples. Not one.

So, from the point of view of visualisation, you could show your first order finite difference as "occuring" between samples ($+\frac{1}{2}Ts$, where $Ts$ is the sampling period). The second order finite difference is something that "occurs" between pairs of samples of the first finite difference. The third order finite difference... (and so on).

Notice here that as you are working "up" the finite differences order, the starting point is "sliding" and the total length of the finite difference signal becomes shorter. The first order finite difference is between $x[n],x[n+1]$, the second order finite difference is between $x'[n],x'[n+1]$ which originate from $(x[n], x[n+1])$, $(x[n+1], x[n+2])$, the third order finite difference is between $x''[n], x''[n+1]$ which originate from $(x[0], x[1])$, $(x[1],x[2])$, $(x[3],x[4])$ and so on.

Notice here that throughout this nested evaluation of these finite differences that "I don't know" bit is becoming larger and larger. You cannot evaluate the $y'''[0]$ of the third finite difference if you don't know the values of the first four samples of $x[n]$. But if you consider an "overlapping" segment across a part of the signal that is long enough to evaluate all finite differences, then you would be able to "see" all of them.

Hope this helps.


Assuming you're looking for symmetry around $x=0$, you can decompose any function into a symmetric and an antisymmetric part:

$$ f(x) = \frac{1}{2}\left(f(x)+f(-x)\right) + \frac{1}{2}\left( f(x) - f(-x) \right)$$

Calling $f_+(x):=\frac{1}{2}\left( f(x) + f(-x) \right)$ and $f_-(x):=\frac{1}{2}\left( f(x) - f(-x) \right)$ you can easily verify that $f_+$ is symmetric while $f_-$ is antisymmetric.

Now you can measure symmetry by taking any suitable norm on your functions. For example the symmetry of $f(x)$ can be measured as

$$S[f]:=\frac{||f_+||}{||f_+||+||f_-||}$$

which ranges from $0$ (fully antisymmetric) to $1$ (fully symmetric).

A possible useful norm on functions is defined as

$$||f|| = \sqrt{\int_\mathbb{R} f(x)^2 dx} $$

which can be discretised to

$$||f|| = \sqrt{ \sum_i f(x_i)^2 \frac{x(i+1)-x(i-1)}{2} }$$


Viterbi Decoder - How does inverting the input bits affect the result?

Convolutional codes are linear in the Galois{2} (i.e. binary) field.  We know this because the convolution operation is linear.  The fact that Viterbi codes are linear has a few implications.  First, if you input all 0's you will get all 0's out.  Second, if you add two inputs together the output will be the sum of their individual outputs.

The last piece of information that we need for a very informal "proof" is that when you input all 1's the output is all 1's (EDIT: this is not always true- see the edit at the bottom of the answer).  We see that this must be so because the encoder state is constant (you are, after all, inputting all 1's) so the output must be constant.  The output cannot be 0's because that would match the all 0's case which would make the code non-linear (other than in the trivial case where the encoder always outputs 0 no matter what the input is).

Now that we've established that inputting all 1's causes the output to be all 1's we can prove that yes, inverting the input is equivalent to inverting the output.

$$
Input => Output\\
Input + Ones => Output + Ones
$$

If we assume the first statement then the second statement follows due to linearity and all 1's producing all 1's, with the understanding that binary inversion is equivalent to adding 1.  0 + 1 = 1, 1 + 1 = 2 = 0.

EDIT: Dilip is correct.  For some generator polynomials an all 1's input will produce an all 1's output and so inverting the input is equivalent to inverting the output.  For some generator polynomials, though, an all 1's input will not produce an all 1's output, so inverting the input is not equivalent to inverting the output.  The generator polynomials that do produce an all 1's output have odd Hamming weights for both polynomials.  The reason that Dilip's polynomial does not meet the criteria is that one of the polynomials ($1 + x^2$) has an even Hamming weight.

EDIT 2: Dilip makes a valid point that even polynomials with odd Hamming weights will produce some zeros at the beginning while the zeros in the state memory are shifted out and replaced with ones.  This is relatively brief, but it does mean that the inversion equivalence does not hold during that time.  It does hold during the "steady state" behavior.


I splitted the question into smaller bits to allow more precise answers.

Question 1: you seem to not understand what spline interpolation is. The quoted sentence means that you need to use an interpolation function that guarantees the continuity of 1st and 2nd order derivatives. This is the case of cubic spline and a justification for their use.
Question 2: yes, I guess that there is a typo here and $n$ is $m$.
Question 3: it is an interval of sampling rates. It is not very clear though...
Questions 4/5: the paper is clearly not clear or informative at this point. What the authors do is to detect a heart cycle by looking at the 1st order differential points (not defined, but that are opposed to the number of local extrema). My guess is that these points are the zero-crossing of the signal derivative (computed from the spline interpolation, cf. question 1).

As an additional note, the paper is not really well written or informative. For example, Laplacian eigenmaps are poorly exposed, and prior knowledge is required to understand it. Maybe you should pick another reference paper in order to get started in this field.


Using Discrete Wavelet Transform in the time-frequency domain

Generally below intensity of the color will tell you about importance of a frequency. Some programming languages have this functionality built in. Let's take for example Mathematica (has free trial).  Let's import some sample recording:

voice = ExampleData[{"Sound", "Apollo11ReturnSafely"}]

Here is built in Spectrogram function:

s = voice[[1, 1, 1]]; r = voice[[1, 2]];
Spectrogram[s, SampleRate -> r]

with time on horizontal and frequencies on vertical axes. For wavelets though it is more natural to use Scalogram. Here is built-in function WaveletScalogram:

cwt = ContinuousWaveletTransform[voice, GaborWavelet[6]];
WaveletScalogram[cwt, ColorFunction -> "AvocadoColors", ColorFunctionScaling -> False]

Scalogram can provide useful info. For example it is great for visualizing voice versus noise feature. Voice is more rich and irregular in structure, noise is more monotonic and repetitive. Horizontal axis is time. Vertical axis is octaves. Higher frequencies are resolved at lower octaves and lower frequencies at higher octaves.

For DWT you will have a similar plot:

cwt = DiscreteWaveletTransform[voice];
WaveletScalogram[cwt, ColorFunction -> "AvocadoColors", ColorFunctionScaling -> False]


Implementing Frequency Domain Convolution in MATLAB (Convolution Theorem)

You need to set the length as well in your fft command.
The two signals are of length $5$ and their convolution is of length $5+5-1=9$. So use this instead:

a=[1 2 3 4 5];
b=[6 7 8 9 10];
x=fft(a,9);
y=fft(b,9);
ifft(x.*y)

ans =

Columns 1 through 8

6.0000   19.0000   40.0000   70.0000  110.0000  114.0000  106.0000   85.0000

Column 9

50.0000

conv(a,b)
ans =

 6    19    40    70   110   114   106    85    50


Extracting top & bottom fields from a YUV 4:2:0 Frame

Assuming the encoder properly low-pass filtered the chroma before sampling (which may not always be true), chroma samples from the adjacent macroblocks might need to be used to properly reconstruct (interpolate) any per pixel color.  Those samples may exist due to motion compensation from previous fields, or might need to be inferred by reflecting, duplicating or extrapolating edge data.


What is the frequency of the Gaussian white noise model in control systems?

No such thing as a single frequency of the noise. That's exactly why it's called white; it has power in all frequency ranges, but not at a single frequency.

Finally, is there a frequency-domain representation of Gaussian white noise?

Yes, a constant power spectral density for all frequencies. That's like white light (which contains also a continuum of all frequencies), which is where it got it name from.
Dealing with continuous spectra, like that of white noise, but also those of actual signals is absolutely bread and butter of control engineering – we very rarely have an input to a system that's always constant or just a sinusoid, which means that its spectrum is never the actual one tone you mention in your question.
I'm not quite sure how to proceed here – I can't really paste three thirds of "Otto Föllinger – Regelungstechnik" (my German standard control engineering textbook) as an answer, and I'm pretty sure that would also not be super helpful to you – because probably, you'll touch upon how you can work with signals and systems in Laplace or Fourier domain, where this really becomes easier.
I'll thus venture a quick "outlook":
Your control system is linear (nod here, please, if you agree, and if you're not certain what linearity is, look it up), thus, if you input a sum of two different signals, the output is equal to the sum of outputs you would have gotten if you used them individually.
Thus, you can start approximating your white noise. Noise has a measurable power. So, pick one arbitrary frequency and put a tone of the same power here. A single line in frequency domain is not a good approximation of a noise that has no lines but is a constant "floor", but at least we've got the same power.
To make the approximation better, you now split the power (i.e. half it, which means multiplying the amplitude with $\frac1{\sqrt2}$) and make two tones that divide the frequency span you're looking at evenly. You calculate the output of your system for both inputs, add it, and get a possible output.
Small problem: You need to assume some phase of your tones. Well, noise can (and will) have independent random phase at different frequencies, so you calculate the output power for a large set of uniformly chosen random phases. You calculate the average power and claim it's an estimate for the output power frequency distribution for white noise.
Instead of halving the power and giving it to two tones, you could of course also had divided the power by 1,000,000 and given it to 1,000,000 tones. Same principle applies. Your estimate of the output power distribution will probably be much better than with 2 strong tones, since 1,000,000 weak tones with random phases are much more alike white noise than 2. However, your math senses are tingling intensely, and you recognize that you're basically doing a large $\sum\limits_f$ over 1,000,000 frequencies and really, you need to do it over uncountably many frequencies, and that sounds like a job for an $\int\cdots\,\mathrm df$.
Thus, you end up with an integral transform, such as the Laplace or Fourier transform. With these, you can describe signals that are deterministic, but have a continuous spectrum, which your Bode plot method seemingly cannot! Think about it: What is the "frequency" of "Input is 0 for eternity, until it turns 6 for a duration of 2, after which it's 0 again for eternity"? There is no single frequency in aperiodic things!
(and by the way, that's why I wrote it has no power at a single frequency: you divide a finite amount of power by an uncountably infinite amount of frequencies. Every frequency gets exactly 0 power. If you integrate over a range of frequencies, even an extremely small one, you get some power, so the result of your division is some kind of power density!)
The fact that the input is stochastic means you're in trouble if you assume fixed phase relationships, and that's why we usually have to resort to power spectral densities. However, from a simple "how does the output power look like if I feed white noise into a low-pass RC filter with a bode plot that looks like ...", you'll see that this is a very useful info to have!
Don't be afraid: The math you've learned for single-frequency inputs is absolutely necessary to build the continuous spectrum understanding, and it will be used all over the place in your continued studies of control engineering (or in fact, all disciplines of electrical engineering, if that's what you're studying). More often than not, the problem really boils down to "well, then just multiply the Bode plot point-wise with the input spectrum and sum / integrate things up"; and that's not really hard to do.


FFT-based coarse carrier recovery for QPSK

A square operation creates an unmodulated tone for a BPSK signal at 2x the carrier frequency (a pure tone for the case that the signal was unfiltered or rectangular pulses with perfect phase and amplitude balance in the BPSK modulation, and typically a stronger carrier with weaker sidebands in the more common filtered or pulse-shaped cases). For QPSK signals (including OQPSK) a fourth law operation is required which produces the same at 4x the carrier frequency.

This can be explained by knowing that multiplying signals causes their phases to add. So a signal multiplied by itself (square) that is modulated 0° to 180° would then become 0° to 360° which is the same as 0° (unmodulated). For QPSK the states are 0°, 90°, 180° and 270° where it can be seen that multiplying any of those phases by 4 results in 0° using the same modulo 360 operation when adding phase.

It is typical when using this approach for carrier recovery to also use a PLL to track the recovered 4x carrier in order to clean up the residual (much lower) sidebands resulting from the finite bandwidth modulated signal and then the PLL-filtered signal is frequency divided to complete the carrier recovery.

I tend not to use this approach given the higher sampling frequency required and the comparative complexity to other methods when implementing all digital systems but this would make for a relatively simple analog approach using analog frequency multipliers. I detail another approach for carrier recovery specific to BPSK, QPSK and QAM that would be in my opinion much more efficient for a digital radio implementation in this post: High modulation index PSK - carrier recovery

This is also an interesting approach by fred harris for a frequency lock loop that I believe would work for QPSK carrier acquisition with additional subsequent phase tracking: How does this FLL work?


Is there two version of ulaw and alaw encoding in audio compression?

Yes. You can choose the µ parameter freely in µ-law encoding. So, there's as many different µ-law encodings as there are possible µ; of course, only very few are standardized.
The same applies to the A in A-law.
You might want to read the wikipedia articles on µ-law and A-law, instead of just the G.711 article.


In general, it means several waveforms will coexist in 5G ecosystem.

In the context of LTE-based 5G (up to now), several subcarrier spacings (hence CP) coexist. The first and foremost motivation for this is to support not only mobile broadband but also other use-cases such as mMTC and URLLC.

Some of benefits of mixed numerology are flexible scheduling and better spectral efficiency.

Anyway, I do agree to other comments that the question is a little too broad.


what are the applications of fourier transform in modulation

OFDM is usually modulated by using an IFFT (or FFT), using the DFT basis vectors for the multi-carrier orthogonality.


How to plot the root of polynomial in Matlab?

If the vector a contains the polynomial coefficients, you can compute the roots of the polynomial by

r = roots(a);

And then you can plot those roots as points in the complex plane with some marker (e.g., 'x'):

plot(r,'x')


How to amplify regions with specific color in an image?

Let us assume you are doing a binary classification with the desired color region set as 1 after applying your algorithm.
Then you can provide a set of color images (histograms) and desired outputs. The weights can be posed as the unknown quantities for this set of linear equations. Lo you have a overdetermined system & you can apply any method like Moore-penrose inversion to estimate the weights.


Why is a circular mask appropriate for Fourier filtering rectangular images?

For simplicity, let's not do any shifting and only consider non-negative frequencies.
Let's assume that the horizontal and vertical image dimensions are even integers $W$ and $H$. Looking at the output of a $H\times W$ 2-d DFT of the image, the $u$th column with $u\le W/2$ represents a horizontal frequency of $u/W$ times the horizontal sampling frequency and the $v$th row with $v\le H/2$ represents a vertical frequency of $v/H$ times the vertical sampling frequency. For a square image grid the horizontal and vertical sampling frequencies are equal and in the following denoted by a single variable $f_s$. The frequency-magnitude of a 2-d frequency at bin $u, v$ will then be $\sqrt{(u/W)^2 + (v/H)^2}f_s$.
For a cut-off frequency $f_c$ your mask would select frequencies:
$$\sqrt{\left(\frac{u}{W}\right)^2 + \left(\frac{v}{H}\right)^2}f_s < f_c\tag{1}$$
$$\Rightarrow\frac{f_s^2}{W^2f_c^2}u^2 + \frac{f_s^2}{H^2f_c^2}v^2 < 1.\tag{2}$$
That indeed defines an ellipse in coordinates $u, v$.
However, if you consider the actual frequencies $\frac{u}{W}f_s$, $\frac{v}{H}f_s$ as coordinates, then what you have is a circle:
$$\text{Eq. 1}$$
$$\Rightarrow\frac{1}{f_c^2}\left(\frac{u}{W}f_s\right)^2 + \frac{1}{f_c^2}\left(\frac{v}{H}f_s\right)^2 < 1.\tag{3}$$
To summarize, it depends on how you express your frequencies.


Partitioning the frequency axis with shifting resolution

After some thinking I managed to create an answer that works for my problem. It does not center around the frequency f0 automatically but there is 4 design parameters that you can tinker with. My plan is to run the design parameters through an optimizer. So I can find the parameters that minimize the difference between the sum of the energy in each partition.

N_partition = 20;                                   % number of partitions
N_switch = -2;                                      % switches how how many partitions there are in
                                                    % the lower vs the upper part of the signal
N_lower = N_partition/2 + N_switch;                 % number of partitions in the lower part of the signal
N_upper = N_partition/2 - N_switch;                 % number of partitions in the upper part of the signal
diff_lower_upper = 0.05;                            % design value that decides how big the first partition in the signal is
a_upper = 25;                                       % design value that decides the slope appearance for the lower part of the signal
a_lower = 2;                                        % design value that decides the slope appearance for the upper part of the signal
vec_u = 0:1/N_upper:1;                              % vector for the lower part
vec_l = 0:1/N_lower:1;                              % vector for the upper part

part_upper = ((a_upper.^vec_u-1)/(a_upper-1));      % the change vector for the upper part of the vector, exponential
part_lower = ((a_lower.^vec_l-1)/(a_lower-1));      % the change vector for the lower part of the vector, exponential
part_lower = fliplr(part_lower).*diff_lower_upper;  % the change vector but corrected for the diff_lower_upper design value

change_vec = [part_lower(1:end-1) part_upper(2:end)]; % the change vector, removes the zero value from part_upper and part_lower

nl = change_vec./sum(change_vec);                   % normalize the change vector, this is our normalized length vector

For my case, this creates the following accumulated frequency/partition. I wanted high "resolution" around my f0 that is 1000Hz and then the partitions should grow in size the farther away they come from f0 which they do. Thanks to everyone that commented!


Fourier Transform of ECG signal in Python

Ok, so a little trick and the transform would be correct. The problem is on last 6th line. Instead of np.fft.fft(original), it should be np.fft.fft(np.array(origitnal).flatten()). I think that the wrong results were because original was a 2D numpy array.


You are already just low-pass filtering the magnitude of the input signal. Right now, your envelope detection "circuit" is just a single-pole filter. Let $f[n]$ be the input to your system and, for simplicity sake, assume that attack_coef $=$ release_coef $ = a_1$. Then we can express the output $y$ of your system as follows:

$$y[n] = a_1 \left(y[n-1] - |f[n]|\ \right) + |f[n]|$$

Which we can rearrange as a transfer function $h$ as follows:
$$y = a_1 \left(y\,z^{-1} - |f|\right) + f$$
$$h = \frac{y}{|f|} = \frac{1 - a_1}{1 - a_1 z^{-1}}$$

You can use this information to either adjust your coefficients to get the desired cutoff frequency, or you could add another pole if you want a steeper roll-off.


The DCT for a signal $x(t)$ on $[0,T]$ is roughly the same as the DFT for the function $x(|t|)$ on $[-T,T]$, that is, mirrored at the vertical axis. Since that function is even, all sine coefficients will be zero and can be removed from the result.

Now $\sin(|t|)$ is different from $\sin(t)=-\sin(-t)$ for negative $t$, so the resulting coefficients will be quite unrelated. If the sampling interval ends at $t=T=2\pi$, the DFT will only exhibit one peak around frequency $f=\frac1{2\pi}$.

But because of the kink at $t=0$ in $\sin(|t|)$, and a second at the end of the sampling interval, the DCT will exhibit a general background of size $f^{-2}$ in the coefficient for frequencies $f>0$.


Is there a function in MATLAB similar to $\tt freqz$ but for $H(\omega)$?

For continuous-time systems and the s-plane where a point $i\omega$ (or $j\omega$ in electronics engineering notation) on the imaginary axis corresponds to frequency $\omega$, there is freqs, see MATLAB documentation.


sketching magnitude of frequency response of H(z)

Since you have the transfer function factored into poles and zeros, it's easiest to draw the pole-zero diagram to get an idea of the overall shape of $H(z)$. The first thing you notice is the pole at $z=1$, i.e. at DC. So your filter is unstable and $|H(e^{j\omega})|$ goes to infinity at DC. Furthermore there's a complex conjugate pole pair on the imaginary axis, i.e. at half the Nyquist frequency. So you can expect a small bump around that frequency. Finally, I would calculate the value of $H(z)$ at Nyquist, i.e. at $z=-1$:

$$H(-1)=\frac{3\cdot 0.5\cdot 0.1}{2\cdot 1.49}\approx 0.05$$

A more detailed sketch than that is indeed difficult, but that's probably not the idea of the exercise.


harmonic waves as integer multiple in spectrum

Every periodic signal can be described as sum of "basis" functions. These basis functions are the sin and cosines at frequencies of integer multiples of a fundamental freqeuency. This fundamental freqeuency is a consequence of the periodicity of the signal. So if a signal is periodic with time $T$, then the fundamental frqeuency is $F_o = \frac{1}{T}$ Hz. That's the basic Fourier series theorem for periodic signals.

It is the waveform of the signal which will define which of these sines and cosines are present.

Consider the figure below:

It is not straightforward to look at the periodic signal above and tell that it will be comprised of so many frequencies. But it is, that's just a consequence that periodic signals are made of a sum of weighted sinusoids at integer multiples of a fundamental frequency.

Just as you require basis vectors in a vector space to define "any signal" in the vector space, the contributions of these basis vectors can vary depending on the vector being defined.Similarly periodic signals live in a vector space whose basis vectors are these sinusoids (that are defined by the periodicity) of the signal and hence there are contributions (big or small) from all these frequencies to define the periodic signal.

So your observations regarding the frquencies appearing in the spectrum are justified. Infact from the fundamental freuency of the motor (150Hz) we even know how fast the motor is rotating :)


Trying to translate formula for blending mode

You can find many ready solution in the Internet. For example here following macro is given:

#define ChannelBlend_Overlay(A,B)    ((uint8)((B < 128) ? (2 * A * B / 255):(255 - 2 * (255 - A) * (255 - B) / 255)))

Additionally you might want to check this package: blImage for source code as you can very easilly translate it. I am sure you will find what you are looking for.


Is ripples in frequency response of a filter ,result of periodic extension of a signal in time domain thereby producing discontinuities?

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


Equalizing a temporal sound signal with a convolution

First, as long as your filter FFT has hermitian symmetry, the convolution with a real signal should also be real.  You may get some residual noise in the imaginary, but you can just pitch it.  If you are getting signal in the imaginary, your implementation is goofed.

Second, what you see is correct.  Because the time signal is real, the frequency signal has symmetry.  The opposite is also true.  Sample zero is the point about which the symmetry occurs.

Third, you have two lengths to take into account.  The sample window, N, and the filter kernel, M.  The sample window will determine how much real time lag you get through the operation, but is inconsequential in rendered applications.  The kernel length determines the number of “taps” in the filter bank.  The padding should be such that each signal is M+N-1 long when the operation is performed.

Fourth, this implementation is flawed because you are directly manipulating the signal in the frequency domain.  What you would want to do is something like this:  every N samples get M samples of the signal.  Perform an FFT.  Do whatever you need to to calculate the filter frequency response.  Set the phase to be linear and centered about the mid point.  Calculate the inverse FFT to get the filter kernel in the time domain.  Covolve the signals (w/ FFT convolution or directly).

Lastly, as I stated in response to your last question, this is not a good implementation of an analog prototype.  Analog filters have poles and should be implemented using biquad filters.  See RBJs biquad filter cookbook for a great source on implementing analog prototypes of common audio filters.  You can still do it the way I’ve described, but to relate it to the analog domain is a misnomer.


Causality, Linearity, and Time Invariance for Systems Described by Linear Constant Coefficient Differential Equations

Note that the system is not uniquely characterized by the linear differential operator. You need auxiliary conditions, and dependent on these conditions, the system may be linear, causal and time-invariant or not. So if you just use the operator $L$ to show that the system is linear then you're actually not showing anything because $L$ is a linear operator by definition.

The output of a system described by a linear constant coefficient differential equation can be split up into two contributions: the zero-state response (ZSR) and the zero-input response (ZIR). The ZSR is the response of the system with zero initial conditions, and, consequently, the ZSR is fully determined by the input signal. The ZIR response is the response of the system caused by non-zero initial conditions, and it is independent of the input signal. Since non-zero initial conditions cause an output term that does not depend on the input signal, scaling or shifting the input signal will not affect that term. Consequently, a system with non-zero initial conditions can neither be linear nor time-invariant, because if it were, the output should scale and shift according to the scaling and shifting of the input signal.

Causality simply means that the system's response at a given time does not depend on the future of the input signal. Clearly, causality is a property of the system, so it cannot depend on the input signal. A causal linear time-invariant (LTI) system is characterized by an impulse response $h(t)$ that satisfies $h(t)=0$ for $t<0$, and it is irrelevant whether the input signal is zero before a certain time or not.

Also take a look at this related answer.


How the zero - phase filter without filtered signal truncation at the end can be implemented?

You can either gather more data at the end than you need for output, or zero pad past the end of your data with at least half the length of your zero-phase filter.

If you can’t get data past the “ends”, and don’t like the artifacts from zero-padding, then you might try other reasonable assumptions (circularity, continuity, etc.) depending on your model.


Isolation of a signal in a frequency band

The algorithm I referenced is premised on the signal being a pure real tone, in which case it will give an exact answer.
The presence of noise or other tones will distort the answer somewhat, but the formulas are actually quite robust.  As long as the tones are fairly far apart (at least several bins), the distortion due to other tones is small.
Aligning your DFT frame so the main tone is roughly a whole number of cycles and a half has two advantages.

The parameters are read from four values of the DFT instead of two, so inherently more noise resistant.

The second harmonic will fall near a whole number of cycles and thus not have much leakage.

If necessary, your signal can be filtered (as you have done) to reduce the magnitude of the harmonics, but you need to know the frequency response of the filter to adjust your amplitude calculation.
Real world signals tend not to be pure tones.  The smaller the segment, the more it will resemble a pure tone.  The formulas I referenced (and am the discoverer of) have the advantage of being exact, so a large N is not necessary for better results.  In a recent job, I was able to take 8 point DFTs on ~1 1/2 cycle frames to calculate the frequency, phase, and amplitude to about 4-5 significant figures.
In radar jobs in the past, and the non-radar one just mentioned, the variation in the signal meant the results you get are dependent on your frame selection and how you choose to average your results.  The algorithm can give you nearly instantaneous values along your signal duration.
The derivation of the frequency formula can be found here:

Two Bin Exact Frequency Formulas for a Pure Real Tone in a DFT

This is an extension to three bins, which can be used if your peak is near a bin (< 0.15 bin width) for more noise robust results.

Improved Three Bin Exact Frequency Formula for a Pure Real Tone in a DFT

In the pure noiseless single tone case, either formula will yield an exact answer.


Why do waveforms that are symmetrical above and below their horizontal centerlines contain no even-numbered harmonics?

The complex Fourier coefficients of a $T$-periodic function $f(t)$ are given by

$$c_n=\frac{1}{T}\int_{0}^Tf(t)e^{-j2\pi nt/T}dt\tag{1}$$

with

$$f(t)=\sum_{n=-\infty}^{\infty}c_ne^{j2\pi nt/T}\tag{2}$$

The coefficients with even indices are

$$\begin{align}c_{2n}&=\frac{1}{T}\int_{0}^Tf(t)e^{-j4\pi nt/T}dt\\&=\frac{1}{T}\left[\int_{0}^{T/2}f(t)e^{-j4\pi nt/T}dt+\int_{T/2}^{T}f(t)e^{-j4\pi nt/T}dt\right]\\&=\frac{1}{T}\left[\int_{0}^{T/2}f(t)e^{-j4\pi nt/T}dt+\int_{0}^{T/2}f(t+T/2)e^{-j4\pi nt/T}e^{-j2\pi n}dt\right]\\&=\frac{1}{T}\int_{0}^{T/2}\left[f(t)+f(t+T/2)\right]e^{-j4\pi nt/T}dt\\&=\frac12\frac{2}{T}\int_{0}^{T/2}\left[f(t)+f(t+T/2)\right]e^{-j2\pi nt/(T/2)}dt\\&=\frac12 d_n\tag{3}\end{align}$$

where $d_n$ are the complex Fourier coefficients of the $T/2$-periodic function $g(t)=f(t)+f(t+T/2)$. These coefficients can only be zero if $g(t)=0$, i.e., if $f(t)=-f(t+T/2)$. This latter condition is exactly the symmetry condition you mention in your question. Consequently, the even Fourier coefficients are zero if (and only if) $f(t)$ satisfies

$$f(t)=-f(t+T/2)\tag{4}$$

Or, in other words, a $T$-periodic function $f(t)$ has only odd harmonics if it is $T/2$-antiperiodic.
I.e., if you shift the function by half its period and flip it across the horizontal axis, it must look the same as before.


Algorithm to identify pitch-shifted recordings?

There is an algorithm called Dynamic Time Warping that allows you to measure a distance between 2 time-series indepentently of their relative time stretching. Thus with a correct threshold value, you could detect if the two recording (the original and the pitch-shifted) are actually the same.

Then, by using the warping path, you could determine the relative amount of pitch shift between the 2 recordings.


How to do the Fourier Transform of bounded function?

Well, imho, one good way to do it is to take the signal's derivative in time and make use of the property $$\frac{d}{dt}x(t) \longleftrightarrow j2\pi f X(f)$$

I'll try to show it analytically.

Since the signal is discontinuous at $t=\pm 2$ you should include Delta function in the derivative, that is, $$\frac{d}{dt}\Big[\frac{A}{2}t\:\mathrm{rect}\Big(\frac{t}{4}\Big)\Big]= \frac{A}{2}\mathrm{rect}\Big(\frac{t}{4}\Big) + \frac{A}{2}t\:\mathrm{rect}'\Big(\frac{t}{4}\Big)$$

Given that a unitary rectangular pulse can be written as a sum of two unit step functions $$\mathrm{rect}\Big(\frac{t}{T}\Big) = u\Big(t+\frac{T}{2}\Big) - u\Big(t-\frac{T}{2}\Big)$$ then the pulse's derivative is $$u'\Big(t+\frac{T}{2}\Big) - u'\Big(t-\frac{T}{2}\Big) = \delta\Big(t+\frac{T}{2}\Big) - \delta\Big(t-\frac{T}{2}\Big)$$

Going back, we get $$\frac{d}{dt}x(t)= \frac{A}{2}\mathrm{rect}\Big(\frac{t}{4}\Big) + \frac{A}{2}t\:\mathrm{rect}'\Big(\frac{t}{4}\Big) = \frac{A}{2}\mathrm{rect}\Big(\frac{t}{4}\Big) + \frac{A}{2}t\Big(\delta(t+2) - \delta(t-2)\Big)$$

The sampling property of a Delta function yields

$$t\delta(t \pm t_0) = t\Big|_{t=\pm t_0}\delta(t \pm t_0) = (\mp t_0)\delta(t\pm t_0)$$

and in your case

$$\frac{A}{2}t\Big(\delta(t+2) - \delta(t-2)\Big)=-A\delta(t + 2)-A\delta(t - 2) $$

and in total

$$\frac{d}{dt}x(t)= \frac{A}{2}\mathrm{rect}\Big(\frac{t}{4}\Big) + \frac{A}{2}t\:\mathrm{rect}'\Big(\frac{t}{4}\Big) = \frac{A}{2}\mathrm{rect}\Big(\frac{t}{4}\Big) -A\delta(t + 2)-A\delta(t - 2)$$

The Fourier Transform of the latter is

$$F\Big\{\frac{d}{dt}x(t)\Big\} = 2A\mathrm{sinc}(4f) - Ae^{j2\pi 2f} - Ae^{-j2\pi 2f} = 2A\mathrm{sinc}(4f) - 2A\cos(2\pi 2f)$$

and from the derivative property $$X(f) = \frac{1}{j2\pi f}\Big(2A\mathrm{sinc}(4f) - 2A\cos(2\pi 2f)\Big) = \frac{jA}{2\pi f}\Big(2\cos(2\pi 2f) - 2\mathrm{sinc}(4f) \Big)$$

which gives

$$X(f) = jA\Big(\frac{2\cos(2\pi 2f)}{2\pi f} - \frac{\sin(4\pi f)}{(2\pi f)^2} \Big)$$

Finally, with the $j\omega$ notation

$$X(j\omega) = jA\Big(\frac{2\omega\cos(2\omega)}{\omega^2} - \frac{\sin(2\omega)}{\omega^2} \Big) =jA\Big(\frac{2\omega\cos(2\omega) - \sin(2\omega)}{\omega^2} \Big)$$


There is not one correct answer to this question, other than the theoretical one which would be that the maximum downsampling factor $M$ is given by

$$M=\frac{f_s}{f_c}$$

where $f_c$ is the highest frequency of the (ideally band-limited) signal. Since there's no ideally band-limited signal, it depends on the application how much aliasing can be tolerated. The frequency $f_c$ must be chosen as the frequency above which the signal (after low pass filtering) has no significant energy. What 'significant' means is determined by the application.


Characterization of the amount of modulation in a signal

As far as the equivalent of an "envelope detector", consider using a moving average of the absolute value (as it would be simple to implement, if you like the approach more elaborate filters could be used). This is how you would do that in Matlab:

M = 10
out = filter(ones(M,1), 1, waveform)

Where M is the number of samples that you want to actually do the moving average over. Note that you will be making a distinction between what you consider the frequency components of your signal vs the frequency components of the envelope as this algorithm seeks to distinguish the two, with the envelope being a low frequency component.

To determine gap duration, you can simply set a threshold on the envelope based on what you decide to be a gap or not and then create a bipolar two level signal (+1 or -1) to measure the duration of the gaps:

gapTest = abs(out- threshold)

As for setting the threshold itself, I would be tempted to use an "N-Sigma" appraoch, where the threshold is determined based on the standard deviation of the data (so would scale automatically and be consistent across sets).

And the following will report the length as "durations" and quantity of gaps as "count" (note that I did it such that leading and trailing gaps are ignored but you may decide to do that differently):

# find edges
x = diff(gapTest)

# ignore leading gaps (as they are partial)
fall = find(x == -1)
x2 = x(fall(1):end)

#find gaps and duration
rise = find(x2 == 1)
fall = find(x2 == -1)

#ignore trailing gaps (as they are partial)
len = min(length(rise), length(fall))
durations = rise(1:len) - fall(1:len)
count = length(durations)

You will need to test the above for all cases and modify accordingly.

As far as the equivalent of a "mod depth" consider using the mean of the data above your threshold compared to the mean of the data below your threshold.


What are the rough steps in obtaining and processing an image from an infrared camera?

It's hard to give a lot of useful information without knowing the details of the signals you are getting from your hardware. Typically the output from any camera is an image, so you would need to decide what is wrong with the image before trying to come up with image processing techniques to solve that problem, you don't seem to have a clear problem specification at the moment.

  What is the difference between preprocessing and processing?

Semantics. It doesn't refer to anything specific in the signal processing domain.

  Is there a preferred order to performing this steps?

Yes but you need to decide which are necessary before worrying about the order.

  Are some absolutely necessary, while others not so much and can be skipped?

It depends entirely on your source data. If your IR camera outputs a clear crisp image that contains all of the information you need then none of the steps will be necessary. If the image is completely washed out then you will need to do contrast enhancement, otherwise it can be skipped. Really nothing useful can be said without more information about your data.

  How would this differ between processing normal light vs IR radiation?

The usual caveat applies that it depends on your data but in principle there is no difference, you are just capturing a different part of the light spectrum.

  Is that roughly the correct chronologically description of the whole process?

No not really, the process will be informed by your source data and your desired output. It's unlikely that all of those topics will apply, more likely one or two of them will and not necessarily in any particular order.

If you want more concrete advice I suggest you explain your procedure in more depth, with the specifics of the data you are acquiring and the results you expect.


It looks like you are plotting a histogram, which is basically an estimate of a probability mass or probability density function. Maybe this is not correct. If your x-axis label is wrong, and this is truly a signal as a function of time, then what you are looking for is an estimation of the actual signal given additive noise. A simpler solution is a frequency-based low-pass filter. People can provide more answers if they know the nature of the signal in your graph.


frequency/phase offset estimation in satellite DOA beamforming applications

Typically, you distribute a mixing oscillator to all ground stations, so that their respective phase is known and can hence be corrected upfront. Having known phase between these receiver parts typically requires a lot of RF design effort and external phase calibration.


How many things can they be done with a spectrogram?

A spectrogram is a reversible (and generally redundant and complex) transformation that is used to separate components, unmix latent signals, some of which being related to frequency analysis (Nota: sometimes, this name is used for the (squared) magnitude of the above. I do use it as a short-hand for the breeds of short-term Fourier transforms).
We can note that the classical Fourier transform is a (simplistic)   form of spectrograms. So anything you can do with Fourier, you can expect to do it (possibly better) with spectrograms.
For instance, spectrograms can be used  to sparsify representations (reduce the proportion of useful components) or to emphasize features from non-stationary signals. So any signal processing technique, that could use the former sparsifying preprocessing (could be an FFT alone)), fits the picture: time-frequency cancellation, detection, noise estimation and removal, segmentation, deconvolution, adaptive filtering, source separation, learning, etc.
For instance, even mp3 or JPEG compression are based on a sort of (non-redundant) spectrograms.


Value of power spectral density $N_0$ or effect of scaling bandwidth to SNR

The noise power continues to be $N_0/2$, independent of bandwidth. The reason is that the noise variance at the output of a filter with frequency response $H(f)$ is $$\sigma_n^2=\int_{-\infty}^\infty \frac{N_0}{2}|H(f)|^2\,df.$$ (This is a direct application of the Wiener-Khinchin theorem). Since you're assuming an orhonormal basis, the noise power is always $N_0/2$, independent of the bandwidth.

Then, doubling your signaling rate without a corresponding increase in spent energy results in halving the SNR.

Regarding your second question, $N_0/2$ is the PSD of the white Gaussian process assumed to exist at the input of your measuring devices. The PSD is indeed measured in W/Hz.

Note that it is a bit weird to define $P$ as watts per second. One watt is one  joule spent per second; watts per second would be joules, or pure energy. Usually, you assume you have a source providing $P$ watts; this source operates for as long as necessary (even from $t=-\infty$ to $t=\infty$).


Do signals with a Fourier transform with discontinuities or zero amplitude (in some frequencies) have Laplace transform?

First of all, it is important to distinguish between the unilateral and the bilateral Laplace transform. For causal signals we can use the unilateral Laplace transform. If for such a signal the Fourier transform exists, then also its Laplace transform exists. This is simply the case because here the Fourier transform is equivalent to the Laplace transform for $s=j\omega$, and if the Fourier integral converges, the Laplace transform must exist with a region of convergence $\Re\{s\}\ge a$, for some real-valued $a\le 0$.

If a two-sided signal is considered, we have to use the bilateral Laplace transform. Here the problem is that for $\Re\{s\}>0$ the signal is damped for $t>0$ but amplified for $t<0$ (and vice versa). Consequently, for two-sided signals with a constant envelop (e.g. sinusoids), the Laplace transform does not exist, even if the Fourier transform exists, because the Laplace integral cannot converge for signals with an exponentially increasing envelope. The same is true for all periodic signals (because they are sums of sinusoids). But even for two-sided signals which decay for increasing $|t|$ the bilateral Laplace transform might not exist if they decay too slowly, which is the case for the sinc function (decaying only as $1/t$).


FIR filter not eliminating out of band noise?

You can solve your problem by an fftshift of the impulse response, as this line shows:

h = fftshift( real( ifft(H)) );

now it should work.

NOTE that there will be a transient of length about $60$ samples (half the linear phase FIR filter length) at the beginning of the convolution output, after which the steady state result emerges...


phenomenon when window size for STFT is too short

It might mean that if you take a window of not enough length so as to contain at least one period of a sinusoid, and instead you are only getting a lot less inside the window of analysis (Like 10% of the period or something like that), then the fourier transform will show you the spectrum of a signal as if that little piece you got inside the window, extended periodically was the entire signal.

To know exactly what is going on in your graph, we need to know the scale of the time axis, and which signal was being analyzed, if it was a sinusoid, what was the frequency. but that might be what the book is trying to say.


Since you're working local it is suggested for you to compare things normalized to their relative surroundings.

The gradient is a vector (2D vector in single channel image).
You can normalize it according to the norm of the gradients surrounding this pixel.

This is what suggested above.


Efficient computation of dot product of convolutions

Convolution is a specific type of linear operator. So in all generality, you can write your expression as $\langle F_1 s_1 | F_2 s_2 \rangle$ with two linear operators $F_{1,2}$. Using the rules of operators in bilinear forms, we can rewrite this as $\langle s_1|F_1^\dagger F_2 |s_2\rangle$ and that again as $\langle s_1|\left(F_1^\dagger F_2\right) s_2\rangle$.

In your case, the two operators have a pleasant property. They are toeplitz, and if we take your padding into account, even circulant. The product of two circulant operators is again circulant, as is the conjugate transpose. That means you can write the result of the operator product as a single convolution. Your convolution kernel is the correlation product $\mathrm{corr}({f_1}^*,f_2)$.

If you implement this using fast convolution, you need to pad this correlation product and the size of your kernel will increase to the sum of the sizes of $f_{1,2}$. Accordingly, the required padding for $s_2$ will also increase.


What is the first derivative of Dirac delta function?

If you imagine a Dirac delta impulse as the limit of a very narrow very high rectangular impulse with unit area centered at $t=0$, then it's clear that its derivative must be a positive impulse at $0^-$ (because that's where the original impulse goes from zero to a very large value), and a negative impulse at $0^+$ (where the impulse goes from a very large value back to zero).
Actually, the generalized derivative of a Dirac delta impulse $\delta(t)$, denoted by $\delta'(t)$, is a generalized function (distribution) with the following properties:
$$\begin{align}\int_{-\infty}^{\infty}\delta'(t)f(t)dt&=-\int_{-\infty}^{\infty}\delta(t)f'(t)dt=-f'(0)\tag{1}\\\delta'(t)f(t)&=f(0)\delta'(t)-f'(0)\delta(t)\tag{2}\\\int_{-\infty}^{\infty}\delta'(\tau)f(t-\tau)d\tau&=f'(t)\qquad\qquad\textrm{(convolution)}\tag{3}\end{align}$$
Property $(1)$ is basically the definition of the derivative of a distribution. The right-most equality is of course a consequence of considering the special distribution $\delta(t)$.
Property $(3)$ means that convolution with the derivative of a Dirac impulse results in the derivative of the convolved function. I.e., the distribution $\delta'(t)$ is the impulse response of an ideal differentiator.
From property $(1)$ (with $f(t)=1$) it follows that
$$\int_{-\infty}^{\infty}\delta'(t)dt=0\tag{4}$$
Furthermore, from $(1)$ it also follows that $\delta'(t)$ is odd, because for every even function $f(t)$ whose derivative exists at $t=0$ we have
$$\int_{-\infty}^{\infty}\delta'(t)f(t)dt=0$$
since for even (and differentiable) $f(t)$, $f'(0)=0$ holds.


Difference between single sided and double sided amplitude spectrum

For baseband real signals the frequency spectrum is conjugate symmetric; i.e., $$X(f) = X(-f)^*$$ This translates to an even magnitude spectrum; i.e., $$|X(f)| = |X(-f)|$$ and an odd phase spectrum.

Because of the symmetry it's also called as the double side band spectrum in commnnications terminology. For the real basband signal, a single side band is also sufficient to represent it. THere are different approaches to obtain the single side band as the uppersideband and the lower side band.

One method of obtaining the upper side band is to convert the signal into the analytic signal which is $x_+(t) = x(t) + j \hat{x}(t)$ where $\hat{x}(t)$ is the Hilbert transform of $x(t)$. The resulting spectrum is:
$$ X_+(f) =
\begin{cases}
2 X(f) ~~, &\text{ for} ~~ f > 0 \\
0 ~~, &\text{ for} ~~ f < 0 \\
\end{cases}
$$

Therefore for the given signal $x(t) = \cos(2\pi 20000 t)$ whose spectrum is $$X(f) = 0.5 \delta(f-20000) + 0.5 \delta(f+20000)$$ the upper side band will be given by $$X_+(f) = \delta(f-20000)$$

Note that you can adjust the amplitude of single side band if you wish to do so.


Does accurate translation estimation using a camera matrix mean the orientation estimation will also be accurate?

If you pose estimation algorithm carefully minimizes the reprojection error using the standard non-linear techniques e.g., Levenberg Mardquardt, and your corner detection scheme is sufficiently accurate in the subpixel level, then yes. I would expect the residuals to be very small. Of course take into account the distortion parameters; they have to be calibrated as well. If you are using OpenCv and referring to the standard functions such as findExtrinsicCameraParam, then you are good to go.

The problem is also known as PnP due to the correspondences. Lepetit et. al. review it nicely here.

Also, Grest et.al provide an overview of 3 methods here.


Polyphase decimation filter plot in Octave

I think I might have finally solved it myself.

Afaik, while it is true that the actual allpass filters are 1st order, the trick is to see that the original formula actually just says that the two chains of allpasses run interleaved. That is the reason for the z^-2 and the reason for the extra multiply by z^-1.
Hence it is perfectly reasonable to simply use z^-1 filters in the implementation, and feed the allpasses the data interleaved.

However, for plotting the response, I had to first compute the cascaded filter coefficients, and then finally apply the z^-1 delay to the A1 chain.

This is the result:

Which is exactly what I wanted.
This is the final code:

c = [0.036681502163648017, 0.13654762463195794, 0.27463175937945444, 0.42313861743656711, 0.56109869787919531, 0.67754004997416184, 0.76974183386322703, 0.83988962484963892, 0.89226081800387902, 0.9315419599631839, 0.96209454837808417, 0.98781637073289585 ];
K = length(c);
i=0;

b0=[1];
a0=[1];
b1=[1];
a1=[1];

while (i<K)

  beta0 = c(i+1);
  b0=conv(b0,[beta0,0,1]);
  a0=conv(a0,[1,0,beta0]);

  beta1= c(i+2);
  b1=conv(b1,[beta1,0,1]);
  a1=conv(a1,[1,0,beta1]);

  b1z1 = [0,b1];

  [A0,w] = freqz(b0,a0,512,2);
  [A1,w] = freqz(b1z1,a1,512,2);

  A=(A0+A1)/2;

  plot(w, 20*log10(abs(A))); % amplitude plot in decibel
  hold on;
  ylabel("Gain (dB)"), axis([0 1 -120 20])
  xlabel("Normalized frequency");
  grid on

  i+=2;
end

I think this should be the accepted answer, but if anyone has anything to add, please do.


What is the effect of sampling time, sampling frequency and SNR on latency?

...clock would vary for each device.

Unless you can afford high end studio devices that can synchronise to a common clock, this is taken for granted. However, saying that, you also need to balance your requirements. If you are using a relatively good "audio interface" (i.e. not a random off the shelf sound card) then you can expect a pretty accurate clock. Beyond this, the question would be: What is it that you do that requires audio synchronisation down to $\frac{1}{48kHz}$ seconds? I am taking 48kHz as a minimum for an audio recording here.

  How will this affect the latency?

The latency of each audio interface remains the same as before.

  Will it increase or decrease? How will the factors like sampling rate, sampling time, SNR and voltage reference impact latency?

For the latency to increase, you would have to introduce some sort of relationship between the audio interfaces. Instead, both are controlled by the same source and as far as each other is concerned, they operate in isolation. All they receive is a command ("start recording") and from there onwards, it is up to the computer/drivers/OS to be fast enough to serve the streams of data produced by each.

Sampling as far as each one of the interfaces is concerned is not affected. For the system overall, you might run into problems if you request the shortest buffer with the highest sampling rate but this then is up to the system the interfaces are attached to. For a given system you might have to try a few different settings to find the better one.

It is interesting that you mention SNR and voltage reference. The SNR of the signal path has nothing to do with the clock path, so, no influence there. Typically, the reference voltage (presumably of an ADC (?)) is inaccessible to you, so, again, no influence there.

If you are so worried that you might get some sort of phasing in the background because of the synchronisation, keep a "beeper" handy (like a clapper, but just for audio and even more accurate). After you hit "record", sound the beeper so that it is recorded through all channels (or at least through two groups that belong to the different interfaces you are trying to sync). Then, when the time comes to put everything together, just make sure that the beeper signals are aligned as best as possible in the beginning of the recording. From that point onwards, assuming good clocks for the individual interfaces, the recordings should remain "locked".

Hope this helps.


$\tt awgn$ and $\tt wgn$ functions in MATLAB always output the same signal, no matter what I input as SNR

So basically you have discrete samples of an audio signal previously recorded and you want to add noise to it. More specifically, you want to add white Gaussian noise (WGN) with a given power to obtain a noisy signal with a specific SNR. That's what the MATLAB function awgn do. You can use it this way:

noiseSignal  = awgn(signal,snr,'measured');

The input snr should be un dBW and the 'measured' config allows MATLAB to automatically measured the power of the input signal (signal should be a vector). The function generated then noise samples (complex or real depending on the signal input), with a power adapted to the SNR value specified, and add them to the signal samples. If you need the noise samples for a reason, you can simply do noise = noiseSignal - signal;.

For more details on the awgn function, you can look at the official documentation.

Furthermore, the noise is generated with pseudo random MATLAB functions, so if you want to always have the same noise samples you will need to have the "same" randomness each time you run the script. For example, you can do that by adding at the beginning of your script,

s = RandStream.create('mt19937ar','seed',529558);
prevStream = RandStream.setGlobalStream(s);

and restoring the default randomness by adding RandStream.setGlobalStream(prevStream); at the end of your script.

Finally, typical values for SNR in dB are integers between [-10;30]. A negative SNR just means that the power of the noise is higher than the power of the signal. The values you took for the SNR (0.004 or 0.9 dB) are really closed, and you will not be able to see the difference on the graphic (but if you look at the noise samples you should). Then, based on your code, changing the power variable is useless since the addnoise function seems to modify the noise power to match the SNR value.

Hope it helped.


How to perform convolution using polyphase structure

I am not quite sure, if i understand you question properly, it would be better if you ask with detail and maybe with a figure, nevertheless if you want convolution of co-efficents of the transfer function then below might be of help. :

and with regard to splitting.. you can easily do it with something like this :

x(1:2:end) = cos(2*pi*f*t(1:2:end));
x(2:2:end) = cos(2*pi*f*t(2:2:end));


Encoding YUV Color Space with Positive Numbers

Why don't you use YUV with a shift?

To make my example simpler, let's use the Color Space YCoCg.
So the channels Co, Cg are in the range [-0.5, 0.5].

When you encode them, make sure you encode Co', Cg' where Co' = Co + 0.5 and Cg' = cG + 0.5.


I have decided to take a different approach: Rather than passing the raw accelerometer data through a filter and choosing between overshoot and ringing versus bias, I am now detecting bias by periodically comparing extrapolated $\Delta \psi _{gyro}$ to $\Delta \psi _{GPS}$:

Initially assume zero bias.
For each extrapolation step, subtract bias from measured $\omega _\psi$
On each new GPS bearing, calculate the change in bias as $\frac{\Delta \psi _{gyro} - \Delta \psi _{GPS}}{t}$
Add this change to the previous bias estimate
Determine the accuracy of the new bias estimate, considering the accuracy of $\psi _{GPS}$, approximate standard deviation of gyroscope output at rest and time elapsed
If we have a previous bias estimate (other than the initial assumption of zero), run the old and new bias estimate through a Kálmán filter; else simply continue with the new bias estimate
Return to step 2.

Results look good so far – this visibly reduces bias without introducing any new artifacts. Responsiveness still needs to be improved, which I can probably tackle by tweaking my accuracy model.


Finding the time domain equation corresponding to a given DFT

(1) I'm not sure my intuition is right

Almost, but you're missing the fact that if you're sampling at 4000kHz, the signals alias; after sampling at 4000kHz, a sine wave at 1500Hz is indistinguishable from one at 2500Hz or 5500Hz, etc.  Also, because the work is being done in complex numbers, negative frequencies have meaning: $e^{-j \omega t}$ is different from $e^{j \omega t}$

  (2) I don't know how to prove it

Go back and study the Fourier Transform more.  You want the Euler identity tattooed onto your frontal lobes: $e^{j \theta} = \cos \theta + j \sin \theta$.  You also want the two offshoots of that: $$\cos \theta = \frac{e^{j \theta} + e^{-j \theta}}{2}$$ and $$\sin \theta = \frac{e^{j \theta} - e^{-j \theta}}{j2}$$

prove these to yourself.  It'll help you out immensely later.

So (after you review the Fourier transform) you have a signal with a real part and an imaginary part: $$\sqrt{32} \left ( e^{j 2 \pi \frac{3}{8}n} + e^{j 2 \pi \frac{5}{8}n} \right) - j\sqrt{32} \left ( e^{j 2 \pi \frac{3}{8}n} - e^{j 2 \pi \frac{5}{8}n} \right)$$  Note that $2 \pi \frac{5}{8} n = 2 \pi (1 - \frac{3}{8}) n$ and recall that $e^{j 2\pi + \theta} = e^{j \theta}$; after some clever trig identities do these look like anything I've just shown you?

  (3) I don't see how to calculate the amplitude A0

It should drop out of the calculations above.

  (4) and finally, is finding x(t) really the answer to the question of writing "the time-domain equation for the discrete () sinusoid in trigonometric form"

Yes, but mind your aliasing!


What are real-valued and complex signals and why is the Fourier transform of a real-valued signal Hermitian?

When one analyses data, outcomes both depends on the nature of the data, and their handling (eg a transformation, like the Fourier transform). One can consider binary data, categorical data, data in two dimensions, etc. To develop a theory and the related tools, it is useful to define where the d ata dwell.

A real signal is simply a signal that takes values that are real. The signal can be discrete or continuous, mono- or mulitdimensional, the observed values are real. The set of real values $\mathbb{R}$ is quite interesting. Since it is a continuous field, it allow a lot of "signal-like" mathematical operations: addition, product, convolution, and a lot of properties related to limits. $\mathbb{R}$ is one of the simplest mathematical structure allowing such operations.

But a natural tool to study signals, spectra,  is related to Fourier analysis. And Fourier uses complex arguments. Why shall we use complex transforms to study real signals? That is a complicated story, but let us say this is mandatory, to cope with linear and time-invariant systems. See for instance: Why cosine is not an eigen signal?

So one can also use complex-valued data. The complex field $\mathbb{C}$ is, maybe, the second simplest mathematical structure allowing the same operations as $\mathbb{R}$, without pain. And a real signal is just a subset of complex-valued data: take a complex object, whose imaginary part of zero, and you have a real signal.

For the last question: real signals are just a subset of complex signals. So when you analyze them (the reals) with a fully complex transformation, you can expect that the outcome will be restricted too. One of the restriction is that real signals have certain symmetries in the Fourier domain, that fully complex signals don't have


How do we know the receiver delay in an OFDM transceiver system?

In basically all simulations I've done, you just send known (but, depending on what you want to simulate, seemingly random) data. Then you correlate the received bits with the transmitted ones, and look for a peak.

That's your end-to-end bit delay, which you don't really care about, but it's how you know which bit you compare to which.

If you data has no significant autocorrelation, and assuming your transceiver basically works with a BER much smaller than 0.5, this is pretty reliable.


If it is difficult for you to remember or calculate the convolution of two sequences then you may try doing it as polynomial multiplication.
Think of x[n] and h[n] as polynomial coefficients. So we have
Px = 3x^2 + 2*x + 1
Ph = 1x^2 - 2*x + 3

Remember that linear convolution of two sequences is polynomial multiplication. Therefore
Py = Px * Ph
Py = (3x^2 + 2*x + 1) * (1x^2 - 2*x + 3)
Py = 3x^4 + 2x^3 - 6x^3 + x^2 + 9x^2 - 4x^2 + 6x - 2x + 3
Py = 3x^4 - 4x^3 +6x^2 + 4x + 3

Now writing the polynomial coefficients back in sequence format
y[n] = {3, -4, 6, 4, 3}


how to apply Short-time fourier transform to code in matlab

this is very simplistic, but it

inputs a wav file and uses only one channel of it
breaks it up into windowed frames,
transforms it to the frequency domain,
does nothing to it in the frequency domain,
transforms back to the time domain,
overlap-adds that frame back into a time-domain array, and
eventually writes it back out to a different .wave file.

try it out.  can't guarantee there ain't typos.

if ~exist('inFile')
    inFile = 'input.wav';
end

if ~exist('outFile')
    outFile = 'output.wav';
end

if ~exist('frameWidth')
    frameWidth = 4096;          % size of FFT frame, better be a power of 2
end
frameHop = frameWidth/2;

analWindow = hanning(frameWidth);

[inBuffer, Fs] = wavread(inFile);

x = [inBuffer(:,1); linspace(0, 0, frameWidth)'];                   % use left channel only, zeropad one frame at the end

clear inBuffer;

numSamples = length(x);

numFrames = floor((numSamples-frameWidth)/frameHop);

y = linspace(0, 0, numSamples)';

n = 0;                              % init sample pointer.  unlike MATLAB, i like counting from 0

for frameIndex = 1:numFrames

     xWindowed = x(n+1:n+frameWidth) .* analWindow;     % get and window the input audio frame

     X = fft(fftshift(xWindowed));              % do the FFT

     Y = X;                         % copy the input spectrum to output

                                % do whatever processing to Y that you like

     yWindowed = fftshift(real(ifft(Y)));           % convert back to time domain, toss the imaginary part

     y(n+1:n+frameWidth) = y(n+1:n+frameWidth) + yWindowed;

     n = n + frameHop;
end

wavwrite(y, Fs, outFile);


Does this lowpass filter have a noncausal impulse response?

Somewhere in the chain between the sound card analog input and the samples that you're plotting, there is certainly a lowpass filter. There is likely an analog anti-aliasing filter before the ADC; in addition, there are likely one or more lowpass filters applied during resampling processes either on the card or in the operating system's audio driver stack (resulting in a stream of samples at the sample rate that you requested).

As you noted, lowpass filters often have an impulse response that resembles a $\text{sinc}$ function. When you excite a linear system with an impulse, the resulting output is just a copy of the system's impulse response (due to the convolution theorem). So, the fact that the observed signal looks like a $\text{sinc}$ function is to be expected. It doesn't sound like this surprised you, either.

However: you did not observe noncausal behavior. If you had, then you should have run immediately to your local patent office. Noncausal systems are not realizable in the real world. Remember the definition: in order for a system's response to be noncausal, its output must lead the input in time. Stated differently, the filter would start outputting its response to the input before you put the input in. Obviously, that's not going to happen.

So what's with the pre-ringing? Simple: the large spike in the middle corresponds to the center of the filter's impulse response. However, that spike does not correspond to $t=0$, the time at which the impulse was inserted into the filter. Consider the following MATLAB example:

% generate a 250th order lowpass filter
b = fir1(250, 0.5);
% plot its impulse response
plot(0:250, b); grid on;

The resulting plot looks like this:

As you can see, the peak in the plot is not at sample index zero. The filter has an overall delay of 125 samples (indeed, all linear-phase FIR filters have a bulk delay of $\frac{N}{2}$ samples, where $N$ is the filter order), and the pre- and post-ringing are centered about that delay. So, when you insert the impulse, you initially see very little response. As the impulse makes its way through the filter taps, the ringing ramps up to a peak at the impulse response's center, then ramps back down to zero.

The takeaway: There is no noncausal behavior in the example you gave. It is possible to simulate noncausal filters in practice by adding enough delay, similar to the delay shown in the lowpass filter's impulse response above.


How to equalize frequency contents using filter?

What you want is not possible, but the problem could be the problem formulation itself. If you have a signal which is non-zero only at a finite number of frequencies, then your signal is a sum of sinusoids (or complex exponentials). What you can do is make their amplitudes equal, but this does not mean that $F(\omega_i)=\alpha$ (where $\omega_i$ are the frequencies of the sinusoids/complex exponentials), because the values $F(\omega_i)$ are not defined. This is because $F(\omega)$ is a sum of Dirac delta impulses, which are only defined by their integrals:

$$F(\omega)=\sum_iw_i\delta(\omega-\omega_i)\tag{1}$$

with some non-zero weights $w_i$. If your goal is to transform (1) into

$$G(\omega)=\sum_i\alpha\delta(\omega-\omega_i)\tag{2}$$

with $\alpha$ being independent of the index $i$, then this is possible by filtering the signal $f(t)$ with a linear time-invariant filter whose frequency response satisfies

$$H(\omega_i)=\alpha/w_i\tag{3}$$

But note that this does not mean that $G(\omega_i)=\alpha$ because the values $G(\omega_i)$ are not defined.


Side Lobes in the Magnitude Response of a Low Pass FIR Filter

The filter is optimized to suppress signal components in the frequency ranges where the side lobes are low. The high side lobes in are just a consequence of optimizing the stopband suppression in those other frequency bands, where aliasing components are expected. The regions with high side lobes are basically "don't care" regions. If you wanted to increase stopband attenuation there, the stopband attenuation in the other - more important - regions would be decreased.


Which kind of machine learning algorithm is Hidden Markov Model?

I would be tempted to reply "none", or "both classification and clustering".

Why "none"? Because HMMs are not in the same bag as support vector machines or k-means.

Support vector machines or k-means are specifically designed to solve a problem (classification in the first case, clustering in the second), and are indeed just an optimization procedure to maximize an "expected goodness of classification" or "goodness of clustering" criterion. The beauty lies in the choice of the criterion and optimization procedure. HMM are not an algorithm per se. They are a specific kind of probability distribution over sequences of vectors - for which we know good parameter estimation and marginal distribution computation algorithms. But asking whether they are in the "clustering" or "classification" family is as ridiculous as asking whether the Gaussian distribution is supervised or unsupervised learning.

Why "both classification and clustering"? Because of the following: Being probability distributions, HMM can be used for classification in a bayesian framework; and being model with hidden states, some latent clustering of the training data can be recovered from their parameters. More precisely:

HMM can be used for classification. This is a straightforward application of the bayesian classification framework, with the HMM being used as the probabilistic model describing your data. For example, you have a large database of utterances of digits ("one", "two", etc) and want to build a system capable of classifying an unknown utterance. For each class in your training data ("one", "two", you estimate the parameters of a HMM model describing the training sequences in this class - and you end up with 10 models. Then, to perform recognition, you compute the 10 likelihood scores (which indicate how likely the sequence you want to recognize has been generated by the model), and the model with the highest score gives you the digit. In the Rabiner tutorial on HMMs, the training stage is "Problem 3", the classification stage is "Problem 2".

HMM can be used in an unsupervised fashion too, to achieve something akin to clustering. Given a sequence, you can train a $k$-state HMM on it, and at the end of the training process, run the Viterbi algorithm on your sequence to get the most likely state associated with each input vector (or just pull this from the $\gamma$ during the training process). This gives you a clustering of your input sequence into $k$ classes, but unlike what you would have obtained by running your data through k-means, your clustering is homogeneous on the time axis. For example, you can extract the color histograms of each frame of a video sequence, run this process on this sequence, and you'll end up with a break-down of the video into homogeneous temporal segments corresponding to scenes (the unpractical bit is that you have to set up the number of scenes $k$ in advance). This technique is commonly used in automatic, unsupervised, structure analysis of video or music.


Calculation frequency response of digital filter with known structure

For the given system you can write down the input-output relation as

$$y[k]=\frac14\left(x[k]+2x[k-1]+x[k-2]\right)\tag{1}$$

because $T$ (or $z^{-1}$) denotes a delay element, which delays its input by one sample interval.
The $\mathcal{Z}$-transform of (1) is (assuming zero initial conditions)

$$Y(z)=\frac14\left(X(z)+2X(z)z^{-1}+X(z)z^{-2}\right)=\frac{X(z)}{4}\left(1+2z^{-1}+z^{-2}\right)\tag{2}$$

From (2) you get the system's transfer function

$$H(z)=\frac{Y(z)}{X(z)}=\frac14\left(1+2z^{-1}+z^{-2}\right)\tag{3}$$

Since the system is stable (any FIR filter is), the frequency response can be obtained by evaluating the transfer function on the unit circle $z=e^{j\omega}$:

$$H(e^{j\omega})=\frac14\left(1+2e^{-j\omega}+e^{-2j\omega}\right)\tag{4}$$

You also could have arrived at (4) by writing down the impulse response from (1)

$$h[k]=\frac14\left(\delta[k]+2\delta[k-1]+\delta[k-2]\right)\tag{5}$$

and taking the (discrete-time) Fourier transform, which also results in (4).

Finally, a few notes concerning misconceptions in your question:

$\mathcal{L}$ usually denotes the Laplace transform, which is defined for continuous functions. $X(z)$ is the $\mathcal{Z}$-transform of the sequence $x[k]$, which can be written as $X(z)=\mathcal{Z}\{x[k]\}$.
The frequency response of a linear time-invariant (LTI) system is the Fourier transform of the system's impulse response: $H(e^{j\omega})=\mathcal{F}\{h[k]\}$, not the Fourier transform of $H(z)$ (which you seem to think, judging from the last equation in your question).


Given a Vertical Derivative Kernel, Find the output image after K is applied to the above input image

Read about convolution here. The link has a web tool to understand convolution better. That will help you. Here is a quick primer on how to go about solving your problem.

You have to choose a pixel, place kernel such that centre of kernel is on the pixel you want to compute. Multiply the overlapping pixels and dive by number of non-zero pixels as shown above. For the boundary pixels you have to take a call as to what you will do. There are many options like use zeros as boundary etc. For derivatives it is better to keep the boundaries as zero. Hope this helped you :)


Why is the Fourier Series a special case of the Fourier Transform and not the other way around?

That's more of a philosophical question. It makes no practical difference if you choose one answer over the other.
Personally, I would agree with the author. I find it  mathematically a lot easier to derive the Fourier Series from the Fourier Transform than the other way around. It also makes more sense to me to say "periodic signals are a subset for all signals" as compared to "all signals are periodic, some just have a period of infinity". That kind of defeats the purpose of having periodicity as a specials case.


Assuming that the system response is smooth, without narrow peaks or narrow valeys. And that the spectrum of the input signal does not vary so much from one chunk to another.
I think this approach is better than doing a single FFT, as you know the FFT of a random signal will have lots of randomness. Averaging may also remove effects due to quantization errors.
Since the FFT bins are equally distributed in the spectrum, each bin corresponds to a very narrow band of your spectrum, and the energy in this band is highly dependent on the signal you used, or noise, if the FFT is big. In order to have a smooth curve you can compute the energy in a wider band by applying a moving average, for instance.


Almost there:

Define a function growRegion(Image, seedX, seedY) that checks the pixel at location seedX, seedY and returns a list of the 4 or 8 pixels surrounding the seed that satisfy the criteria for inclusion. In your case, that is just one threshold so you could easily say, "if the pixel is adjacent to the seed pixel AND it is above the threshold then include it to the list of potential points".
Establish a list of points to be checked. For instance seedPoints=[(seedX,seedY)].
As long as seedPoints is NOT EMPTY

Pop a point (someSeedX, someSeedY) from seedPoints (remove it from the list)
Call the function returnPoints = growRegion(Image, someSeedX, someSeedY)
For each point in returnPoints mark a white pixel on some Mask image.
Extend seedPoints with the result of growRegion.

At this point, your Mask image contains all 4-8 connected pixels that satisfy the criteria for inlcudes. seedPoints fills up 4 (or 8) times the rate at which it empties so for large regions this might seem to run slowly.

Hope this helps.


About using Walsh matrix as spreading sequence

Yes you can, because you do spreading for each bit of data. So, regardless the length of your data, you are going to spread every bit on your code. In case if you have $x = 10$ you will have data after spreading with length 8, and if $x=101101$ you will have data  to send with length 4 x 6 = 24. Just try to understand the process of spreading (I think it's by XOR), then you do spread for each bit of data, not for all data.

Good luck


The minimum level time (i.e. minimum of L and H) must be greater than the sampling period. This way, there is no possibility for a pulse to be missed.

So, $T_s < \rm{min}(L,H)$.


You're right, the first Fourier transform correspondence in your reference is wrong. It should be

$$\mathcal{F}\{e^{-\alpha n}u[n]\}=\frac{1}{1-e^{-\alpha}e^{-j\omega}}\tag{1}$$

You just need to substitute $a=e^{-\alpha}$ and use the formulas you know.


Matlab FFT not producing symmetric spectrum

The problem comes from the DFT size. Try the code below that uses power of 2 fft.

T = 1e-4;
fs = 1/T;
pulse = rcosdesign(.5, 6, fs);
L = length(pulse);
Nfft = 2^nextpow2(L);
ft = fft(pulse,Nfft);
ft = fftshift(ft);
ft_mag = abs(ft);
freq=(-Nfft/2:Nfft/2-1)/Nfft*fs;
plot(freq,ft_mag)
axis([-3 3 0 max(ft_mag)])


Difference between filtering using Matlab and filtering using C code generated by Matlab

You have differences in the transient, but the steady-state looks the same therefore it is probably a memory initialization problem in the C-code

Here are a few things to check :

Are you sure your input is exactly the same?  Make sure your input buffer is properly initialized in the C-code.

You probably have some kind of data buffer to store the previous samples. Again make sure the data is properly initialized. If the data buffer is not properly initialized you will likely have garbage in the data buffer instead of having '0' in the data buffer which would explain why the transient is different.

The Matlab coder generates 3 user functions, Initialize, Step and Terminate. Make sure you call Initialize once, before make any call to Step.


Can we band-limit a PCM signal? Is it even relevant to relate a totally Analog term(Band limit) with digital one(PCM)

Yes -- it's just a matter of selecting the right pulse shape.

(By the way, it makes perfect sense to talk about the bandwidth of a PCM signal; it is, after all, a continuous-time physical signal with bandwidth, energy, power, etc, just like any other signal.)

Most introductory textbooks (irritatingly) assume that a PCM signal must use rectangular pulses. You don't have to, and actually you almost never do, except in the simplest applications.

The PCM signal is a train (or sequence) of pulses with varying amplitudes. This train can be written as $$s_{\text{PCM}}(t) = \sum_k a_k p(t-kT_p),$$ where $T_p$ is the pulse rate, and $a_k$ is the amplitude of each pulse.

In the simplest PCM signal, $a_k$ are either 0 or 1, and $p(t)$ is a rectangular pulse of duration $T_p$. The bandwidth of $s_{\text{PCM}}(t)$ is the same as the bandwidth of $p(t)$, which in this case is (mathematically) infinite.

By selecting a different pulse $p(t)$, you can limit the bandwidth of $s_{\text{PCM}}(t)$. A popular pulse shape is the "raised cosine", which has a bandwidth between $1/T_p$ and $1/2T_p$.

Note that the bandwidth can also be limited by "brute force": just filter the rectangular $s_{\text{PCM}}(t)$ with a low-pass filter of the desired cutoff frequency. This will distort the signal, but in simple applications (high SNR, no further distortion or attenuation) this usually causes no significant performance degradation.


The most naive and easiest method would be to interpolate the data points. If you know what is the maximum frequency present in your signal then you can specify the sampling period. Let's say you pick the lowest possible: 12ms. Now you can perform some interpolation of your data. Normally I would suggest you to perform the sinc interpolation, but you can start with linear/spline/cubic.

In MATLAB/Octave you can use the interp1 function to interpolate the data points. If your data is called x, known time stamps are t, and new time stamps are t_new, then you can do the following to calculate new dataset X:

t_new = 0:0.012:max(t); %spaced by 12ms up to maximum value of vector t
X = interp1(t, x, t_new, 'spline');

This is very simple approach but should work - although you might get some aliasing. You can also check out this article for more intelligent methods: Nonuniform Sampling and Reconstruction in Shift-Invariant Spaces.

Now when you have your dataset, you can perform any filtering you like. Can I ask you what your sensor is so it's providing such data?


Is it correct to assume that $(-1)^n = \cos(\pi n)$ while computing the DTFT of $(-1)^n$?

$$(-1)^n = \cos(n\pi)=(e^{j\pi})^n=(e^{-j\pi})^n=e^{j\pi n} = e^{-j\pi n}$$

All of these would work. However, applying a DFT is bit tricky since it's directly at the Nyquist frequency and it technically violates the sampling theorem. You can certainly do the math but interpreting the result is not straightforward. For example if you phase shift $\cos(n\pi)$ to $\sin(n\pi)$ the magnitude of the FFT drops to zero, although you only time-shifted the signal.


Transposed FIR filter structure to compensate for extra delays on addition

Ah, its easier than I thought, just mirror the introduced delays. So this means introduce two additional $z^{-2}$ in the last two stages of the transpose FIR filter structure. After some basic math to get $H(z)$ the coefficients are as below in bd2, which gives the same magnitude frequency response as the original.

bd2 = [0 0 0 0 b(1) b(2) b(3) b(4)];


PSNR of two images of different size in matlab

When interpolating, why not use the same size as the image to compare with:

% some random inputs
original = rand(435,782);
img = rand(100,200);
% position where to interpolate
[x,y] = meshgrid(linspace(1,size(img,2),size(original,2)),linspace(1,size(img,1),size(original,1)));
% interpolate
dest = interp2(img,x,y,'bicubic');
% display size of result
disp(size(dest));


Adjust Matlab Gui parameter while playing an audio

As far as I know, the simplest, very old, matlab function to play wave PCM audio is "sound(x,Fs,Nbits)" and it is a high level function that will buffer the entire piece before it is played. And it will not allow adjustment on "x" before it is fully played. If your play function is also similar in mechanism, then you shall consider either using a low level audio playback function or divide your audio into small pieces and send each piece in sequence. However this last method may create problems due to multimedia timing issues. And as a result played audio may have pops & cracks in it.