Keeping track of units for spectral energy of discrete signals

First it should be noted that the nomenclature using power and energy as properties of signals is somewhat misleading. For physical energy to come out of a time integral over a squared signal, the signal would need to be of dimension $\sqrt{[\mathrm{energy}]/[\mathrm{time}]}$, which is not the case for practically all signals that you encounter.

In fact, if you have signals that fully describe the state of a physical system over time, and the energy of the system can there be expressed using just that signal, the result is usually very different and most importantly, it is a property of each time instant and not the whole signal. If you have a signal that describes a non-temporal aspect of a physical signal, like the momentary state of a string, the energy of the system will in general require the time derivative of that state.

Therefore, the terms energy and power are a bit of stretch and should be used carefully. I personally prefer speaking of a squared magnitude and a squared magnitude density.

With this preamble, the use of physical units like time and frequency in the definition of power and energy of a signal become somewhat relative. The choice of scaling and therefore of units is mostly arbitrary and depends on what you want to express and which transformations you want to keep your measures invariant.

For example, if you have a choice of sampling your physical system at different rates and you want the energy measure to be invariant under such rate changes, then you need to first understand how the amplitudes behave under a rate change. The normal behaviour of an amplitude is to be be unaffected by a change of sampling rate. Some acquisition methods may however rely on integration to record the samples and your sample values may therefore scale with the sampling rate. So your energy sum may contain different powers of the time step in order to keep the energy invariant under sampling rate changes.

Keeping track of all this for different kinds of signal behaviours quickly become tricky. That is why I made a habit of defining a discrete signal not just by giving a sequence of numbers, but also by providing a homomorphism to an infinite dimensional function space that allows me to work out all the meaningful measures and properties in a well defined framework with properly defined integral transforms and continuous variables.

So, a discrete signal is a tuple $\{ s\in \mathbb{S},\, D: \mathbb{C}\to\mathbb{S},\, \bar{D}:\mathbb{S}\to\mathbb{C} \}$ where $\mathbb{S}$ is the space the discrete signal lives in, $\mathbb{C}$ the space of the continuous time signals, and $s$ a representation of the discrete time signal in $\mathbb{S}$. Note that $D\circ\bar{D}=\mathrm{id}_\mathbb{S}$ but not necessarily $\bar{D}\circ D=\mathrm{id}_\mathbb{C}$.

You give these two maps $D$ and $\bar{D}$ all the invariance properties you would like to see. Afterwards you can pull back all definitions from $\mathbb{C}$ to $\mathbb{S}$ by applying the two domain changing maps.

This method of working with discrete signals will give you a well defined way of talking about properties define for continuous signals and allow you to get your scaling correct under sample rate changes or changes of other parameters relevant for your discrete signal representation.

I will not try to replicate your derivation using this approach, because it is not clear to me which invariance and signal properties you want to preserve and what exactly your goal is. Like I said, it depends on the context what is right. I still hope this answer helps your understanding and answers the slightly more fundamental question of how to deal with units, dimensions and scaling of discrete signals compared to continuous signals.


Applications for the Maximum Phase Filter

we can deconvolve linear-phase filters into a minimum-phase filter and its reverse a maximum phase filter

We can formulate this more broadly. Any LTI system can be split into a cascade of it's minimum phase filter and an all pass (which is indeed a maximum phase filter). So,
$$H(z) = H_{min}(z) \cdot A(z)$$

where $H_m(z)$ is the minimum phase filter that matches the magnitude $|H(z)|$ and $A(z)$ is an allpass , i.e.

$$A(z) = k \cdot \frac{\prod_{n=0}^{N-1}z-1/p_n^*}{\prod_{n=0}^{N-1}z-p_n}$$

The proof for this is straight forward: Let's assume that H(z) has a zero $q_k$ outdside the unit circle. We can simply factor this out as

$$H(z) = H_k(z) \cdot (z - q_k) = H_k(z) \cdot (z - q_k) \cdot \frac{z - 1/q_k^*}{z - 1/q_k^*}$$

$$ = H_k(z) \cdot (z - 1/q_k^*) \cdot \frac{z - q_k}{z - 1/q_k^*} = H_{k,min}(z) \cdot A_k(z)$$

where $H_{k,min}(z)$ is the version of $H(z)$ where the zero $q_k$ has been made minimum phase by inverting it to $1/q_k^*$  with $|1/q_k^*|  < 1 $ and $A_k(z)$ is a first order allpass with a zero at $q_k$. Note that the magnitude stays the same, i.e. $ |H_{k,min}(z)| = |H(z)| $

That type of splitting an LTI system into it's maximum and minimum phase part is can be quite useful.

Inverting: The minimum phase part is invertible and the inverse is also minimum phase. The inverse of the maximum phase is simply it's own time reverse. That's not directly invertible (since the inverse is non causal), but you easily "eyeball" how much bulk delay you need to add to get enough of the phase response correct.

Audio: Linear phase FIR are mostly useless in audio. Getting decent resolution at low frequencies requires many 1000s of taps and the computational complexity and the latency are often prohibitive. Equalizer or cross over design is often split into two steps: minimum phase filter to get the amplitude right and than allpass filters to do time/phase alignment as needed.

Building blocks: Allpass filters are very useful building blocks to build more complicated structures. For example you an build odd order Butterworth cross overs as sum + difference of two allpasses and as a bonus you get a signal that's phase matched to the perfect reconstruction path. Warped FIR filters are another useful topology where allpasses are quite handy. Beamformers need more or less independent control of amplitude and phase.

So in general in an application where you want to control phase and  amplitude more or less independently, splitting into minimum (for amplitude) and maximum (for phase) can be quite useful.


Random process $X(t)$ with autocorrelation function given find the mean and the variance

The limit $\lim_{\tau\to\infty} R_x(\tau)$, if it exists, equals $E^2[X(t)]$ and so $E[X(t)]=0$ in this case.

More generally, the mean of a WSS process is nonzero only if the power spectral density has an impulse at the origin. This can be applied to
periodic autocorrelation functions such as $\cos(\omega_0t)$ pointed out
in @MattL's comment. If the Fourier series for a periodic autocorrelation
function has a nonzero DC term, the mean is nonzero.


In the context of sampling, an image is a copy of the spectrum of the original signal occurring at multiples of the sampling frequency. If the images overlap, then aliasing occurs, i.e. frequency components above $f_s/2$ cannot be distinguished from components below $f_s/2$ (where $f_s$ is the sampling frequency), and the original signal cannot be completely recovered from its samples.


aliasing folding phenomena in frequency domain matlab

Your plot is correct. You are sampling three waves $f_1 = 10\textrm{Hz}, f_2 = 30\textrm{Hz}$, and $f_3 = 70\textrm{Hz}$, with a sample frequency of $F_s = 1.5\times70\textrm{Hz} = 105\textrm{Hz}$. This means that your Nyquist frequency is $F_N = F_s/2 = 52.5\textrm{Hz}$, and corresponds to the maximum value of your frequency axis.

As such, the signals $f_1$ and $f_2$ will be correctly sampled without problems, but the signal $f_3$ is higher than the Nyquist and so should appear as an alias (folded back into the sampled frequency range) at the frequency $f_{\textrm{alias}} = F_N-(f_3-F_N) = F_s-f_3 = 35\textrm{Hz}$. This is what you obtained.

If you want to create a plot similar to the one in the figure, then you could try to first sample your waves with a much higher sampling rate, in order to fully capture them all faithfully. Then, on the same plot, display what happens if you undersample with $F_N<f_3$.


Why doesn't law of large numbers apply to this stationary time-series?

I believe that you are thinking that each value of $X_t$ is determined by a different realisation of $Y$, which in this example is not true.

Suppose that $Y$ is the value that comes out from a dice throw. Thus, $Y$ can achieve any integer value between $1$ and $6$. Throw the dice once. Suppose that you get $2$. Then, according to the definition in the example, your time-series would be:

$$X_t=2 \qquad {\text{for all }}t$$

What I think the quote from Wikipedia is saying is that the mean value of this stochastic process doesn't coincide with the mean value of the random variable $Y$, as

$$\mathbb{E}[X_t]=2 \neq 3.5 = \mathbb{E}[Y]$$

If you constructed the time-series letting each value of $X_t$ be determined by a new throw of the dice, then the means would match.


Compute SQNR (Signal to Quantization Noise Ratio)

The signal voltage is a sinusoid signal variant from $-V$ to $+V$, thus the RMS value

$$V_{rms} = \frac{V}{\sqrt{2}} = \frac{2^N \Delta V / 2}{\sqrt{2}}$$

because $2V = 2^N\Delta V$ where $\Delta V$ is quantization step.

Quantization noise is modeled as uniform random variable in $[-\Delta V /2, +\Delta V /2]$ thus its standard deviation is $V_n = \Delta V / \sqrt{12}$.

SQRN is defined as $\mathrm{SQNR}=\frac{V_{rms}^2}{V_n^2}$

$$\sqrt{\mathrm{SQNR}}=\frac{V_{rms}}{V_n} = 2^N \sqrt{\frac{3}{2}}$$

$$\mathrm{SQNR}_{dB}=20\log_{10}\left(\frac{V_{rms}}{V_n}\right) = 20\log_{10}\left(2^N \sqrt{\frac{3}{2}}\right) = 6.02N + 1.76$$

You must have mistyped your question.


Additive synthesis vs. "just draw the waveform"

If your goal is to get a sawtooth wave, then additive synthesis is a bad way to go about it as a sawtooth wave is a superposition of an infinite number  of sine waves (in the continuous case).  Many instruments, particularly woodwinds, are well modeled by a relatively few sine waves.  The exact mixture determines the timbre and thus what instrument it sounds like, e.g. an oboe as compared to a clarinet.  Violins, however, generate a slip-stick motion that is more like a sawtooth wave and thus are harder to model accurately and efficiently with additive synthesis.

Unless you are going for an "electronic" sound, I would say using LinTable by itself would be a somewhat unusual thing.  However, using it as input to a filter would make more sense.  For example, if you did use a sawtooth wave to model a violin, you'd almost certainly pass it through a (among other things) low-pass filter to round it out a little and change the timbre.


You can try using a Pitch Track for monophonic signal, you need track low pitches!

There is a wide range of algorithms to find Pitches in time domain (zero crossing, AMDF, YIN, MPM, ASDF, autocorrelation ....) or Frequency domain (shs, hps, cepstrum ....)

Each one of them possess advantages or disadvantages !

With the Pitch Track you can convert one frequency to note.

I never used Flash AS3 to DSP, one fast and easy way are the autocorrelation based algorithms.


Pixel Wise Write into a Video File in MATLAB

I assume your matrix is 3D where the first 2 dimensions are Width and Height and the third is time (Gray Scale Video).

If you can write your processing as a Filter you can use MATLAB's filter() function.
The function has the dim property hence you can set it to 3 and it will apply on the filter over the third dimension.
The limitation is applying the same filter to whole data.
If you need to have different filter for each pixel then you should iterate over the filter bank.

Another solution is to use the functions permute() and reshape() to rearrange the data to have the Time Axis in the dimension which fits your needs.


Periodicity property of DFT in time and frequency domain

If the expression that defines the DFT is evaluated for all integers
  $k$ instead of just for $k = 0, \dots, N-1$ , then the resulting
  infinite sequence is a periodic extension of the DFT, periodic with
  period $N$.

  The periodicity can be shown directly from the definition:

  $$ X_{k+N} \stackrel{\mathrm{def}}{=} \ \sum_{n=0}^{N-1} x_n e^{-\frac{2\pi i}{N} (k+N) n} = \sum_{n=0}^{N-1} x_n e^{-\frac{2\pi i}{N} k n} \underbrace{e^{-2 \pi i n}}_{1} = \sum_{n=0}^{N-1} x_n e^{-\frac{2\pi i}{N} k n} = X_k.$$

  Similarly, it can be shown that the IDFT formula leads to a periodic
  extension.

  Source: Wikipedia.

While the periodic property of DFT isn't very widely utilized, it often causes aliasing problems.

Source: http://www.dspguide.com/ch10/3.htm


Mathematical model of equivalent-time sampling, the resulting unevenly spaced periodic signal, and its interpretation

Since the physical measurement takes a certain time $t_s$ (leading to the original sampling rate), each sample in $s[t]$ can be seen as a convolution of a sampling function (perhaps a boxcar) $f(T)$ with the continuous physical value function $s(T)$ during the measuring time. I suspect this is different in the questions that are linked above, which regard economical data and might be regarded instantaneous values.

If that is correct, then in equivalent time the highest frequency that is sampled should not be higher than in the real / original time. Unless there is some way to deconvolve the measurements in equivalent time, the sampling rate in both timeframes should be equal. So the high number of samples between the start and the end of the average period should be reduced to match the true sampling rate, before comparing two measurements.

I think that a moving average with subsampling can be performed. This should be set such that the resulting sampling rate equals the original one. Since in equivalent time the samples are unevenly spaced this can not be performed by a filter in the frequency domain. The moving average has to be performed in the temporal domain instead.

Scale phase $\varphi$  to the average period $P$ of all $N$ periods so its unit is equivalent to that of $t$

$t'(\varphi) = \frac{\varphi}{N} \sum_{i=1}^{N} t_{i,end} - t_{i,start}$ to give $s[t']$

Moving average of $s[t']$ for every $nt_s$

$s'[nt_s] = \frac{\sum_{t'\geq (n-1)t_s}^{t'\lt nt_s} s[t']}{\sum_{t'\geq (n-1)t_s}^{t'\lt n t_s} \frac{t'}{t'}}$

Im not sure about the notation since every $t_s$ has a different number of samples and the period length is not equal to an integer number $n$ of $t_s$. In practice you could append $s[t'+P] = s[t']$ and then use $t_s' = \frac{P}{ceil(\frac{P}{t_s})}$, leading to a slightly interpolated result, but at least it is truly periodic.

The resulting waveform has no more detail / higher frequency content than the original signal. This signal can be compared between experiments like any other signal, so at least it is a workable solution.

This might be a poorly constructed method compared to other suggested solutions (Kalman filter, non uniform FFT´s gridding and unevenly spaced time series operators), but its probably easier too. Not sure if I lose any information in the process.


what is proper setting for Fb,Fc in Complex morlet wavelet (cmor)?

Remember that Wavelet Transforms are nothing but time-localized filtering/correlation operations. The wavelet transforms provide a unified framework for getting around the Heisenberg Uncertainly Principle that the Fourier Transform suffers from. So when you ask "what should my settings be for bandwidth, and center frequency", you are asking for filter parameters to be given to you. No one can give you filter parameters but yourself, and your application. This is why they are parameters. :-)

More concretely, how to pick your $F_c$ and $F_b$ have to do with the nature of your signal. First, some pre-liminaries:

Wavelet Transforms are nothing but matched filters that you construct, such that you have the best match to the 'features' in your signal of interest. One set of 'features' that you might consider would be the neighbourhood of center frequency, and the neighbourhood of bandwidth. Again, no one can pick those for you, because you need to know what you are looking for. However, you can pick a range commensurate with what you expect. (This is actually one of features that make Wavelet Transforms very powerful).

The morlet wavelet is given by the following function:

$$
\Psi(t,f_c,T_p) = \frac{1}{\sqrt{\pi T_p}} \ e^{\frac{-t^2}{T_p}} \ e^{-j2 \pi f_ct}
$$

Here, $t$ is the time index in seconds, $f_c$ is your center frequency, and $T_p$ represents what I have called a "Period Parameter". Inspect the above equation. Notice how the morlet wavelet is in fact nothing but a complex exponential centered at frequency $f_c$, windowed by a zero-mean gaussian function, with $\sigma = \sqrt{\frac{T_p}{2}}$

Thus, to answer your first question, the $F_c$ parameter, must be chosen such that it is at or near the frequency you wish to interrogate.

Regarding the wavelet bandwidth: What I have called $T_p$ here is your "Period Parameter", and this is directly related to your bandwidth. Go back to the Morlet Wavelet equation, and consider a scenario where instead of a gaussian window, we used a rectangular window. In the Fourier Domain, we attain a sinc function due to this type of window, and we can say that the Null-To-Null Bandwidth, $B$, is simply: $B = \frac{2}{T_{period}}$, (in the passband).

The figure below illustrates this for a complex exponential located at $f_c = 2$ Hz, with varying $T_{period}$ values:

You can clearly see that as the period extent is increased, the bandwidth decreases. This is the inverse nature of the time-frequency relationship. If you wanted to use this pseudo-wavelet, what would you pick the $T_{period}$ to be? The answer would be "Whatever $T_{period}$ my signal-of-interest hangs around at". In fact, this is exactly what is done in the Short-Time-Fourier-Transform, (STFT), although one would pick a smarter window instead of the boxcar.

Now let us re-inspect the Morlet Wavelet Equation, and re-insert the gaussian windowing function, with $\sigma = \sqrt{\frac{T_p}{2}}$. If we do this, we get the following time-domain plots, and the corresponding Discrete Fourier Transform (DFT)s of the Morlet Wavelets:

Once again, notice how as the extent - as the variance of the gaussian window shading the complex exponential increases, the bandwidth in the Fourier Domain decreases. (In this case, we have much nicer side-lobe behavior, but this is tangential). In any case, You can once again see how we have an inverse relationship between the $T_p$ parameter, (or the standard deviation parameter of your gaussian, whichever suits you), and the bandwidth of the corresponding filter.

What to do with this knowledge: Now you understand how bandwidth is related to the gaussian shading function embedded within the Morlet Wavelet. Thus, for any signal $x[n]$ you have, you may simply take the DFT of your signal, and observe its bandwidths at various center frequencies. Then, adjust the $\sigma$ and $f_c$ parameters of your Morlet Wavelet accordingly, such that the DFT function of the Morlet filter directly overlaps the DFT function of your template signal. This is how you can set your parameters.

The code to generate the above plots is here:

clear all;
fS = 500;
tStart = -4;
tStop = 4;
timeVector = linspace(tStart,tStop, (tStop-tStart)*fS );
fCenter = 2;
tP = [0.2 1 3];
exps = length(tP);

figure(1); clf(1);
figure(2); clf(2);
for ii = 1:exps

    tPeriod = tP(ii);
    timeMask = zeros(1,length(timeVector));
    timeMask((timeVector >= -tPeriod/2) & (timeVector <= tPeriod/2)) = 1;

    psiFakeWavelet = exp(2*1i*pi*fCenter.*timeVector).*timeMask;
    psiWavelet = ((pi*tP(ii))^(-0.5)).*exp(2*1i*pi*fCenter.*timeVector).*exp(-timeVector.^2/tP(ii));

    %Demonstrating the Rectangular Window
    figure(1);
    subplot(3,2,2*ii-1);
    plot(timeVector,real(psiFakeWavelet), timeVector,imag(psiFakeWavelet));
    ylim([-1.2 1.2])
    xlabel('Time / Seconds');
    title(sprintf('Morlet Wavelet, T_{period} = %2.2f, f_c = %2.2f Hz', tPeriod, fCenter));

    input = psiFakeWavelet;
    Nfft = 10 * 2^nextpow2(length(input));
    psd = 20.*log10(fftshift(abs(fft(input,Nfft))));
    freqs = [0:Nfft - 1].*(fS/Nfft);
    freqs(freqs >= fS/2) = freqs(freqs >= fS/2) - fS;
    freqs = fftshift(freqs);
    figure(1);
    subplot(exps,2,2*ii);
    plot(freqs, psd);
    xlim([-10  15]);
    xlabel('Frequency / Hz');
    title (sprintf('PSD, Null-Null-BW = %2.2f Hz', 2/tPeriod));

    %Demonstrating the Morlet Wavelet
    figure(2);
    subplot(3,2,2*ii-1);
    plot(timeVector,real(psiWavelet), timeVector,imag(psiWavelet));
    ylim([-1.2 1.2])
    xlabel('Time / Seconds');
    title(sprintf('Morlet Wavelet, T_{period} = %2.2f, f_c = %2.2f Hz', tPeriod, fCenter));

    input = psiWavelet;
    Nfft = 10 * 2^nextpow2(length(input));
    psd = 20.*log10(fftshift(abs(fft(input,Nfft))));
    freqs = [0:Nfft - 1].*(fS/Nfft);
    freqs(freqs >= fS/2) = freqs(freqs >= fS/2) - fS;
    freqs = fftshift(freqs);
    figure(2);
    subplot(exps,2,2*ii);
    plot(freqs, psd);
    xlim([-10  15]);
    xlabel('Frequency / Hz');
    title ('Wavelet PSD');

end


The condition in order to not lose information after upsampling and downsampling

The answer to this is to design a filter that perfectly passes the original spectrum and completely rejects the images that the decimation and interpolation process creates. Such a filter is not achievable, so the outcome becomes a design constraint on how much distortion we can tolerate (as in all filter design). Helpful for this is to understand the mechanism that creates the images in this case so I am including some graphics below that I have that further explain this process.

Interpolation by inserting zeros

We upsample with minimum distortion by inserting zeros in between each sample in the time domain, and then follow this with an interpolation filter. The interpolation filter passes the original spectrum and rejects the images the zero-insert creates.

The graphic below shows the image replication of the spectrum by upsampling by 4 (insert 3 zeros between each sample). Note how we started off with just our desired spectrum in the "first Nyquist zone" which is the freqency from DC to half the sampling rate ($F_s/2$) for a real signal. (The same applies to the spectrum for a complex signal, just in that case the spectrum is unique from $F_s/2$ to $+F_s/2$ so we would include that or equivalently the signal from DC to $F_s$.) After the zero insert we end up with 4 additional spectral images in the new first Nyquist zone at the higher sampling rate, for the case of a 4x sampling rate increase. (So in your case you would have N such images).

The ideal interpolation filter would "grow" the zeros to the actual value needed for perfect interpolation. Such a filter could be a multiband filter (Matlab firls and firpm do these easily) as depicted in the graphic below, but a low pass filter would also be fine to use. The actual distortion depends on how much of the image frequencies we allow through and how much pass band distortion (either loss or ripple) we allow, this all becomes a design constraint on the complexity of the filter, tolerable delay etc...

Decimation by dropping samples

For decimation we reduce the sampling rate by simply removing samples (So if we decimate by M, we only use every Mth sample). The spectral challenge of images in the decimation process is very similar. While in interpolation we were concerned about clean-up in that we filter after we do the zero-insert; for decimation we are concerned with filtering prior to the dropping of samples due to aliasing; signal energy that will fold into our spectrum of interest if we don't filter prior to removing samples.  The graphic below shows this for an example of decimating by 4. This is identical to the aliasing issues in A/D conversion (in fact A/D conversion is effectively decimating from an infinite sampling rate! So here we are doing a digital to digital conversion):

Like the interpolation filter, the ideal decimation filter would pass our spectrum of interest with no distortion while completely eliminating any signal that may exist in the higher spectral locations that can fold in. Note if you know by design that these spectral locations have already been filtered (including amplified noise floor), then filtering may not be necessary.

Like the interpolation filter design, a multiband filter approach would provide more rejection for the same number of taps (more efficient):

And finally a rational rate conversion (both decimate and interpolate) suggests that you can do the interpolation first so that the filtering can be shared between the two operations. This is a consideration dependent on the actual filtering design requirements for both operations, but feasible that one low pass filter design could cover both requirements. In some cases with high numbers for M and N it may also makes more sense to break up the interpolation and decimation operations further to simplify the filtering requirements.

In all cases a significant limit to what can be achieved is how much of the digital spectrum available we are using. The more excess spectrum available, the simpler the filter design becomes since the design is driven by how tight the transition bands need to be.


The idea of a spectogram is to split your signal into a number of blocks or frames, which are potentially overlapping. After windowing, an FFT is calculated per frame. The output of these FFTs are collected as column vectors in your graph. Thus, the x-axis is related to time and the y-axis is related to frequency. Since the FFT of a real-valued signal is symmetric, only half the frequencies have to be plotted. Therefore, the y-axis goes from 0 to Fs/2=128.

Since you have chosen the FFT length (NFFT=256) equal to the lenght of your signal, you have only one full frame of data.

Potential solutions are to reduce your FFT length NFFT or increase your data length. For now, the easiest way seems to be in reducing NFFT. Note however, that you also need to adjust noverlap in that case. A widely used amount of overlap is NFFT/2.

Example:
The following command will give you a spectogram.

specgram(data, NFFT=64, Fs=256, noverlap=32)

Alternative implementation:

Nfft=64
specgram(data, NFFT=Nfft, Fs=256, noverlap=Nfft/2)

The FFT length parameter NFFT leads to a tradeoff between time and frequency resolution.

EDIT To answer additional questions in the comments.

From the documentation in
[http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.specgram]:

  Returns the tuple (spectrum, freqs, t, im):

  spectrum: 2-D array

  columns are the periodograms of successive segments

  freqs: 1-D array

  The frequencies corresponding to the rows in spectrum

  t: 1-D array

  The times corresponding to midpoints of segments (i.e the columns in >spectrum)

  im: instance of class AxesImage

  The image created by imshow containing the spectrogram

Thus, the first element in the returned tuple contains the power values.


Identify whether to have unique output in this ARMA system

As I mentioned in my comment, I think there is a typographical error in the equation. I think it is supposed to be
\begin{equation}
y[k] = \sum_{i=1}^{M}a_iy[k-i] + \sum_{j=1}^{N}b_jx[k-j].
\end{equation}
This is an example of a linear constant-coefficient difference equation. These are introduced very early in most textbooks about digital signal processing. For example, they are introduced in Section 2.5 (page 33) in the 1989 edition of Discrete-Time Signal Processing by Oppenheim and Schafer.
Given an input signal $x$, there will not be a unique solution signal $y$. This is because a specific solution requires specifying the input signal as well as initial conditions for $y$. Different initial conditions lead to different solutions.
You may gain some insight into the subject from Oppenheim himself in lecture 3 of a DSP course recorded in 1975 and now posted online.


Sampling Frequency and Spectral Regrowth

Yes that's called as the aperture effect and indicates that if the FFT window length (the aperture) coincides with an integer number of signal periods, then the resulting FFT spectrum will have no leakage artefacts (mostly having zeros, up to quantization noise, bins) and significantly nonzero bins only at the signal frequencies.


Since convolution describes the operation of a linear time-invariant (LTI) system, the question is if the effect of an LTI system can be compensated by another LTI system. In the discrete-time domain you can use the $\mathcal{Z}$-transform to analyze LTI systems. If a signal $x(n)$ (with $\mathcal{Z}$-transform $X(z)$) is filtered by a system with impulse response $h(n)$ (with $\mathcal{Z}$-transform $H(z)$), then the $\mathcal{Z}$-transform of the resulting output signal $y(n)$ is given by

$$Y(z)=X(z)H(z)$$

Recovering $X(z)$ from $Y(z)$ can theoretically be done by

$$X(z)=\frac{Y(z)}{H(z)}$$

except at points $z$ where $H(z)=0$. The above equation corresponds to convolution with an impulse response $g(n)$ with $\mathcal{Z}$-transform

$$G(z)=\frac{1}{H(z)}\tag{1}$$

The question is if this system can be realized by a causal and stable filter. This is indeed the case if $H(z)$ is a minimum-phase system. Assuming rational transfer functions (which we always get when discrete-time systems consist of adders, multipliers, and delay elements), this means that $H(z)$ is a causal and stable system with all its zeros inside the unit circle of the complex plane (its poles are of course also inside the unit circle because of stability). Due the inversion in (1), the zeros of $H(z)$ are the poles of $G(z)$ and vice versa, i.e. if the zeros of $H(z)$ are inside the unit circle, then $G(z)$ is stable. One more thing that can be seen from (1) is that if $H(z)$ is an FIR system, then $G(z)$ is IIR. Nevertheless, in practice equalization if often done by FIR filters, which approximate the inverse filter in some optimal sense.

Let's look at a simple example. Assume an FIR system

$$h(n)=\delta(n)-a\delta(n-1),\quad |a|<1$$

Then we have

$$H(z)=1-az^{-1}$$

$H(z)$ is a minimum-phase system because it has one pole at $z=0$ and a zero at $z=a$, which is inside the unit circle if $|a|<1$.
The inverse system is given by

$$G(z)=\frac{1}{1-az^{-1}}=\sum_{n=0}^{\infty}a^nz^{-n}$$

the impulse response of which is

$$g(n)=a^n,\quad n\ge 0$$

which is of course IIR.

So in this case the convolution

$$y(n)=(x*h)(n)=x(n)-ax(n-1)$$

can be perfectly compensated for by the convolution

$$x(n)=(y*g)(n)=\sum_{k=0}^{\infty}g(k)y(n-k)=\sum_{k=0}^{\infty}a^ky(n-k)$$

Note that exact deconvolution of a minimum-phase FIR system can only be achieved by an IIR system. The only systems that can be inverted (exactly) by FIR systems are all-pole IIR systems, i.e. systems with all their zeros at the origin of the complex plane (like $g(n)$ in the example above). Note however that in practice exact compensation is usually not the goal. E.g., in digital communications, almost all (linear) equalizers are implemented by transversal (FIR) filters, because they can be adapted more easily and because they are always stable.


Mean of innovation sequence is zero if the Kalman filter is optimal. A non zero mean implies that there are certain errors present in your assumptions of covariance of noise or even your system matrices(if you are using a state space model).


Spectral averaging without phase information

The short and most common answer is that you can average without phase but there are cases where you can average phase as well.

Cross correlation between
2 sigmals where one has a delayed component of the other that is of a "long" duration can exploit complex averaging.

Spectral analysis using the periodogram averages magnitude or magnitude squared.

There are numerous examples of where you do or don't use complex averaging.  It depends on the problem being solved and is typically going to be relatively obvious if you should or shouldn't.  If you expect the phase to evolve in a predictable way over your observations, you can use that to your advantage.


Stochastic Methods for Image Deconvolution Problem

Any Deconvolution method which takes into account the noise in the image is basically a stochastic approach.
Usually, the model for Deconvolution is:

So having the noise in there makes it a problem with stochastic properties.
Remark
If by stochastic you meant sampling from the Posterior Distribution then you may have a look at Stochastic Image Denoising by Sampling from the Posterior Distribution (Though it is not about Deconvolution but on Denoising).


How to change the angle unit labels on polar plot in MATLAB?

When you opens the polar plot in the figure editor you can find out that the coordinates around the plot are not from the axis.
It are hidden text objects (Press on one you see it). So this means you have to edit them by yourself. a work around for this is this:

figure
t = 0:.01:2*pi;
polar(t,sin(2*t).*cos(2*t),'--r')

hHiddenText = findall(gca,'type','text');
Angles = 0 : 30 : 330;
hObjToDelete = zeros( length(Angles)-4, 1 );
k = 0;
for ang = Angles
   hObj = findall(hHiddenText,'string',num2str(ang));
   switch ang
   case 0
      set(hObj,'string','0');
   case 30
      set(hObj,'string','-30');
   case 60
      set(hObj,'string','-60');
   case 90
      set(hObj,'string','-90');
   case 120
      set(hObj,'string','-120');
   case 150
      set(hObj,'string','-150');
   case 180
      set(hObj,'string','-180');
   case 210
      set(hObj,'string','-210');
   case 240
      set(hObj,'string','-240');
   case 270
      set(hObj,'string','-270');
   case 300
      set(hObj,'string','-300');
   case 330
      set(hObj,'string','-330');
   otherwise
      k = k + 1;
      hObjToDelete(k) = hObj;
   end
end
delete( hObjToDelete(hObjToDelete~=0) );

If you want other directions etc i think you do that manually otherwise the plot start tot mess up.


The Matlab code of the S-transform can be found here for example: http://en.verysource.com/code/1180181_1/st.m.html

The main advantage compared to CWT is its simplicity. In fact, it can be seen as equivalent to a Morlet wavelet, see for example paper by Ventosa, Schimmel, Simon, Dañobeitia and Mànuel "The S-transform from a wavelet point of view", IEEE Trans on Signal Processing, 2008.


ROC of Z transform Doesnt include a pole on the boundary?

The ROC $0.5<|z|<\infty$ does not imply that there's a pole at $z=0.5$. What it does say is that there's at least one pole satisfying $|z|=0.5$ (and no other pole with a radius larger than $0.5$), and this is also the case for the given denominator polynomial. The roots of the polynomial are
$$z_{1,2}=0.5\, e^{j\pi/3}\tag{1}$$
and
$$z_{3,4}=0.5\, e^{-j\pi/3}\tag{1}$$
I.e., there are two double poles, and all four poles lie on a circle with radius $0.5$.


"Noise Spectral Density" a.k.a "Power Spectral Density" in datasheets

Power spectral density is given in power per $Hz$. "Input Noise Spectral Density" is given in units of volts per $\sqrt{Hz}$.  These are related:

$$
PSD = \frac{W}{Hz}
$$

$$
W = \frac{V^{2}}{R}
$$

$$
Input\ noise\ spectral\ density\ = \frac{V}{\sqrt{Hz}}
$$

The $\sqrt{Hz}$ is from the conversion power to volts.  For example if you want the noise ($nV$) in $1kHz$ $BW$ you must multiply $\frac{nV}{\sqrt{Hz}}\sqrt{1000}$


I'll try to give you some hints to get you started. First of all, you should know that the bilinear transform is given by

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

If the analog prototype filter is normalized such that its cut-off frequency is $\Omega_c=1$, then the constant $k$ in (1) is given by

$$k=\frac{1}{\tan\left(\frac{\omega_0}{2}\right)}\tag{2}$$

where $\omega_0$ is the desired cut-off frequency of the discrete-time filter. Some more background on this is given in this answer.

You should also know that the poles of the normalized analog filter lie on a circle with radius $1$ centered at $s=0$ (of course the poles only lie on the left half of the circle). If the filter order is odd as in your example (why?), there must be a pole at $s=-1$. Now you can try to figure out how that pole is transformed to the $z$-plane, i.e.

$$\frac{1}{1+s}{\huge|}_{s=k\frac{z-1}{z+1}}=\ldots\tag{3}$$

From (3), knowing that that pole is transformed to a pole at $z=0.6$, you can determine the constant $k$, and from that, via Eq. (2), you can compute $\omega_0$. Now you know that the filter's attenuation is $3\text{ dB}$ at $\omega_0$, and it is $44\text{ dB}$ at $\omega=1.2$.

The last thing you should know is that with a Butterworth lowpass filter of order $N$, you get approximately $6N\text{ dB}$ attenuation per octave, at least in the frequency range considered here. Now figure out how many octaves there are between $\omega_0$ and $\omega=1.2$, figure out how many dB's difference in attenuation there must be between those two frequencies, and from this compute an estimate of the filter order $N$.


I hope you do not mind, but I am going to change terminology a little bit. Since you mentioned simulations, suppose you are sampling the Gaussian white noise, at the input to the integrator, at a constant rate of $f_s$ samples per second. The point spacing, $Δt$, is $1/f_s$ seconds. Suppose $N$ independent consecutive samples are collected. Then the $N$ samples span a period of $τ_a$ seconds and this equals $N Δt$. The sample mean of the $N$ noisy inputs to the integrator obviously provides an estimate of the known $V_{in}$ input voltage. The sample variance of the $N$ noisy inputs to the integrator provides an estimate of the known $σ^2$ variance of the Gaussian white noise. The power SNR at the input is $(V_{in}/σ)^2$. The empirical estimate would be the square of the sample mean divided by the sample variance and bigger $N$ values are almost surely better.

The Gaussian white noise has zero mean and $σ^2$ variance, so the bilateral noise power spectral density, $η$, is $σ^2 Δt$. The ideal integrator has $τ_i$ integration time constant. (For an ideal op amp integrator, with just an input resistor, $R$, and feedback capacitor, $C$, $τ_i = RC$. The sign inversion is trivial and can be ignored.) Now suppose that the integration period, i.e., integration aperture or gate duration, is $τ_a$ seconds, exactly as above. Assume the integrator has initial condition zero, integration starts at $t = 0$ and ends at $t = τ_a$, so $τ_a$ is simply the reciprocal of the OP’s $F_s$. N.B. My $f_s$ is not the OP’s $F_s$.

Then the DC power gain of the integrator is $(τ_a/τ_i)^2$ and the noise equivalent bandwidth is $1/2τ_a$. Now consider two cases.

Case 1: Let $τ_i = τ_a$.  Then the signal power at the output of the integrator, at $t = τ_a$, is simply $(V_{in})^2$, because the DC power gain is unity. The noise power is $2η$ times unity DC power gain times $1/2τ_a$, which equals $η/τ_a$. But $η = σ^2 Δt$ and $τ_a = N Δt$. So the output power SNR is $(V_{in})^2 N/σ^2$. This is $N$ times larger than the input power SNR.

Case 2: Let $τ_i ≠ τ_a$.  Then the signal power at the output of the integrator, at $t = τ_a$, is $(V_{in} τ_a/τ_i)^2$, because the DC power gain is $(τ_a/τ_i)^2$. The noise power is $2η$ times $(τ_a/τ_i)^2$ times $1/2τ_a$, which equals $η τ_a /(τ_i)^2$. But $η = σ^2 Δt$ and $τ_a = N Δt$. So the output power SNR is $(V_{in} τ_a/τ_i)^2/[σ^2 Δt τ_a /(τ_i)^2]$, which again equals $(V_{in})^2 N/σ^2$. This is again $N$ times larger than the input power SNR.

In Case 1, the signal power is independent of $N$, while the noise power is inversely proportional to $N$, via $τ_a = N Δt$. So the power SNR is proportional to $N$. This is simply averaging independently sampled white noise. In Case 2, signal power is proportional to $N^2$, while noise power is proportional to $N$. So the power SNR is again proportional to $N$.

Relationship to matched filtering

The figure below shows that, if nothing was known about the input waveform for $t > τ_a$, then the input could be treated as though it was a rectangular pulse starting at $t = 0$ and ending at $t = τ_a$. The integration would effectively be gated integration, which would be a simple case of matched filtering: the impulse response of the (aperture synchronized) gated integrator is a rectangular pulse of duration $τ_a$ and reversing it changes nothing because the gated integrator's impulse response is symmetric. The figure below shows the calculation of the matched filter power SNR. The integration time constant, $τ_i$, has been set equal to unity because it affects both the pulse energy, and the noise, the same way, so it divides out.

I hope this helps!


The notation to express a datum filtered by arbitray number of one-dimensional filters?

...my question is, is there an elegant notation to express such operation?

If you know that these operations are to act on orthogonal dimensions, you can define their 1D version and then use indexing to express the order by which the function is applied to the multidimensional signal.

But if you are working with multidimensional signals over an orthogonal vector basis throughout, you can simply specify this in the beginning and use (for example) $\mathcal{F}(x)$ throughout, to denote a multidimensional discrete Fourier transform.

Hope this helps (?)


Convolution of one signal with an evenly spaced signal

This is known as polyphase decomposition. It is often used as en efficient implementation of filtering combined with decimation or interpolation.


Is there a way to calculate modulation depth after the modulation process?

If you look at the envelope of an AM signal then for 100% modulation the amplitude of the envelope will vary between 0 and 2 * Ac, where Ac is the amplitude of the carrier. It's a linear relationship, so for e.g. 50% modulation the amplitude will vary between 0.5 * Ac and 1.5 * Ac. So as long as you know that you have a full scale modulation signal then you can derive Ac from the max and min envelope amplitude (mean of Amax and Amin) and from there you can calculate the modulation depth (Amax / Ac).


From another stackexchange question (https://ham.stackexchange.com/a/2018) you should use a filter approximately 4 times as wide in Hz as the expected words per minute of the morse code.  That answer also includes some links and more info.

The value "numtaps" says how long to make the filter.  You are using an FIR filter.  The sharpness and the stop band attenuation (how much the stuff you don't want to hear gets reduced) are both dependent (amongst other things) on the value of numtaps.  Higher values gets you better filter performance, but it is more work for the computer.
If you find that you need high values, then you are not using an optimal filter.  There is another factor that influences the quality of the filter, and that is the "window" that is used to shape the filter.  If you don't use a window, then the stop band will not be reduced far enough.

The length of the filter (numtaps) is also relevant to the lowest frequency you can work with.  If you try to make a band pass filter, for example, then numtaps needs to be at least as long as the lowest frequency (in samples.)  Say you are sampling at 1000Hz and want to make a band pass from 1Hz to 3Hz.  You will need for numtaps to be a least 1000.  If you used instead only 333 for numtaps then you would end up with a lowpass that lets through everything below 3 Hz.

You seem to be using scipy - here is the description of the FIR-filter designer "firwin." http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.signal.firwin.html

Also refer to the windows page: http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.signal.get_window.html#scipy.signal.get_window

I think you should try a blackmanharris or some other window that doesn't require parameters - otherwise you will have to start guessing numbers that you don't understand.

Simply using a window at all (other than boxcar) should improve things considerably so that you don't need such long filters - numtaps equivalent to the beep frequency should not be needed - you should be able to get away with well under 100 for numtaps


Filtering of complex exponential functions

HINT:

If the given $X(\omega)$ is the Fourier transform of $x(t)$, what does $x(t)$ look like? It will turn out that removing exponential terms with certain values $\tau_i$ is equivalent to suppressing certain time intervals in the time domain.


Changing limit of integral make time-variant to time-Invariant system

The system-1 whose input/output equation is
$$ y(t) = \mathcal{T}\{x(t)\} = \int_{-\infty}^{t} x(\tau) d\tau $$ is time-invariant.

Whereas the system-2
$$ y(t) = \mathcal{T}\{x(t)\} = \int_{-\infty}^{3t} x(\tau) d\tau $$ is time-varying.

On the other hand the system-3
$$ y(t) = \mathcal{T}\{x(t)\} = \int_{-\infty}^{3t-3t_0} x(\tau) d\tau $$ is still time-varying for any (fixed) value of $t_0$.

I presume that your teacher have tried to underline the fact that what makes the system-2 time-varying is the fact that $y(t-t_0) \neq \mathcal{T}\{x(t-t_0)\}$ which becomes as the equation shows:

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

and
$$\mathcal{T}\{x(t-t_0)\} = \int_{-\infty}^{3t} x(\tau-t_0) d\tau = \int_{-\infty}^{3t-t_0} x(\tau) d\tau $$

These two would be equal if th upper limit of the second equation would be $3t-3 t_0$. That's what your teacher told. Since they are not equal therefore the system-2 is time-varying.

Note that artifically making the limit as $3t-3t_0$ is quite meaningless, as the fixed value of $t_0$ has nothing to do with the variable amount of shifts that occurs in $x(t-d)$ or $y(t-d)$.


What is the meaning of half window functions?

Because the head related impulse response (HRIR) is zero for negative time, you may as well interpret it so that they have used a full Hann window.


Which of these signals have same Nyquist rate as $x(t)$?

Your understanding is right.

With the usual Fourier transform pair notation :

$$x(t) \longleftrightarrow X(\omega)$$

and assuming bandwidth of $x(t)$ is $W$, then the four cases will be

1-) $x^2(t) \longleftrightarrow X(\omega) \star X(\omega) \implies Bandwidth = 2W$
2-) $x(2t) \longleftrightarrow \frac{1}{2}X(\omega/2) ~~~~~ ~\implies Bandwidth = 2W$
3-) $ x'(t) \longleftrightarrow j\omega X(\omega) ~~~~~ ~~~~~\implies Bandwidth = W$
4-) $ x(t)\star x(t) \longleftrightarrow X^2(\omega) ~~~\implies Bandwidth = W$

You can see that, cases 3 and 4 has the same bandwidth as $x(t)$ and therefore their Nyquist sampling frequencies will also be the same.


determine the impulse response of given second order difference equation

The impulse response of a system is its zero-state response to an impulse at the input. If the system is linear and time-invariant (LTI), then the system's response to any input signal can be described in terms of the impulse response. The impulse response only depends on the structure of the system, and it cannot depend on any initial conditions. In fact, a system with non-zero initial conditions is not linear because a part of the output (the zero-input response) does not depend on the input signal but solely on the initial conditions and the system itself.


New to DSP and having trouble filtering data

You should try with implementing your own linear filter. It is not really difficult, and you will have more control about what's going on behind the scene and it will make it easier for you to tune the filter bandwidth.
Speaking of bandwidth, you can try to plot the frequency analysis of your data using Fourier transform, that will help you in finding a suitable cutoff frequency.
A linear filter implements the equations $\sum_{i=1}^n a_i y[i] = \sum_{j=1}^m b_i x[i]$. As long as you normalize the weights $(a_i), (b_i)$ to sum to 1, you should preserve the range of the data.


Why do we care if a system is BIBO stable?

The acronym BIBO, when used by university DSP professors, is merely a fancy (academic) term that means stable. So when when you encounter the term BIBO in the literature of DSP you can replace that term with the word "stable".

As for your final question, the answer is yes. Think of this: You have an audio signal that you'd like to apply to the input of a lowpass filter. And it's possible to design an IIR lowpass filter (a recursive filter, with feedback) using commercial software using filter coefficients that are represented by super-high precision floating point binary numbers. Then when you quantize your coefficients to 8-bit values in order to implement your lowpass filter using an 8-bit microcontroller, it's possible for the "imprecise" 8-bit filter coefficients to cause the IIR filter to become unstable.

That scenario, described in every DSP textbook, means each input sample applied to your filter will cause the filter's successive output samples to uncontrollably grow in value. And when fixed-bit-width binary numbers grow in value they will eventually experience overflow error.  The bottom line here is: When an information-carrying discrete signal is applied to an unstable digital filter the output quickly becomes random noise. To quote my young granddaughter, "This is a bad thing."


Finding Repetetion rate of ns pulsed waveforms

problem statement: my sampling rate would be 2MSPS - to acquire a 100ns pulse i require that much

Nope. 100ns is 5 times as short as the interval between two sampling instances in your sampling device. If it does indeed measure the voltage accumulated over a capacitance, this might work – if that capacitance is as small as possible, which it often is, you might just miss pulses.

However, you're right, the signal you're interested in actually doesn't have a bandwidth of 10 MHz – the maximum frequency of interest is 50 kHz, so observing a system that indicates the presence of a pulse (not the pulse itself) with a 100 kS/s sampling rate should suffice.

If you know you just got presence information (i.e "there was an impulse during the last sample time, or not"), an FFT is a bit over the top. Simply count the samples between an observation.

The hard part actually would be the electronics to your photo diode – 1ns impulses imply that your system needs to have a bandwidth of 1 GHz – so you'd be deep in the world of RF/microwave circuitry. What you'll probably end up doing is amplifying the current coming from your diode, with something avalanche-based, or a sensitive microwave transistor. Don't forget to lay out your board following the needs of a microwave signal! You can't afford spurious interference here.

You'd essentially amplify the energy of the pulse. The voltage over a simple RC lowpass fed with that energy with a 90% discharge time of let's say 1/(100 kHz) would serve as a feasible indicator of activity; observe that, if you want with 2MS/s, and you'd have captured a signal from which you can directly see the impulses.

Notice that photodiodes with a > 1GHz bandwidths will be noisy, so you'll probably have some false impulses. Here, your FFT approach might actually tell the periodicities from the white noise.


Can I create a transformation matrix from rotation/translation vectors?

I was stuck on the assumption that the destination points in the call to getPerspectiveTransform had to be the corners of the output image. Once it dawned on me that the destination points could be somewhere within the output image I had my answer.

float boardX = 1240;
float boardY = 1570;
float boardWidth = 1730;
float boardHeight = 1400;

vector<Point2f> destinationCorners;
destinationCorners(Point2f(boardX+boardWidth, boardY));
destinationCorners(Point2f(boardX+boardWidth, boardY+boardHeight));
destinationCorners(Point2f(boardX, boardY+boardHeight));
destinationCorners(Point2f(boardX, boardY));

Mat h = getPerspectiveTransform(detectedCorners, destinationCorners);

Mat bigImage(image.size() * 3, image.type(), Scalar(0, 50, 50));

warpPerspective(image, bigImage, h, bigImage.size());

This fixed the perspective of the board and everything in its plane. (The waviness of the board is due to the fact that the paper wasn't lying flat in the original photo.)


Recovering an image using IDCT through DCT coefficients

If you used dct2  on an image, or on each image channel separately for color images,  with the appropriate input format, then idct2 recover the original image (the transformation is invertible and orthogonal), up to rounding errors (similar discussion in Why is $\exp(\ln(x))−x\neq 0$ in floating point arithmetic?).

You can find for instance a demo below:

dataRGB = double(imread('autumn.tif'));
data = rgb2gray(dataRGB/255);
dataDct=  dct2(data);
dataDctIdct = idct2(dataDct);

subplot(2,2,1)
imagesc(data);
xlabel('Original (grayscale)');
subplot(2,2,2)
imagesc(log(abs(dataDct)));
xlabel('Original DCT');
subplot(2,2,3)
imagesc(dataDctIdct);
xlabel('Original recovered');
subplot(2,2,4)
imagesc(data-dataDctIdct);
xlabel('Difference');
colorbar
colormap gray


What Does It Mean When Cross Correlation Between Two Images Occur at Negative Coordinates

It means the best match to template happens outside the image.
For instance, let's say your template is 5 by 5.
And you got answer which is -1, -1. It means the part of the image which best matches you image is centered at [-1, -1] and you only have part of it in your image.
This is really an extreme case.
P. S.
If you share your data (2 Images) we'll be able to assist more.


Advantage of applying a Window Function in Analysis and Synthesis of STFT?

The STFT is an instance of an Analysis/Synthesis system. If we restrict to windows of finite length, it amounts to windowing blocks of length $L$, taking the DFT, and hopping by $K$ samples to perform the same operation. An A/S system is said invertible, or perfect, if you can recover the original signal from the blocks of DFT coefficients.
If $K>L$, some samples are missed, there are holes between blocks, and generally some information is lost (unless you have additional information).
If $K=L$, blocks are contiguous, and do not overlap.  Here, the transformation is critical, or non-redundant: you have the same numbers of Fourier coefficients than of samples. In that case, you can invert it even with non-rectangular windows. It suffices that the windows don't vanish, because windows without zero-values can be inverted (and Fourier as well).  There is only one inverse in that case. This is for your first question.
If you allow  overlap $K<L$, then several things may happen. I shall come back later, but basically:

a proper synthesis window can make a system with a vanishing analysis window invertible;
it can reduce artifacts resulting from  performing processing after the analysis stage (frequency removal, thresholding etc.);

In the continuous time domain, if you have an analyzing window $h(t)$, for any synthesis window such that
$$ \int_{-\infty}^{+\infty}g(t)h^*(t)\mathrm{d}t = 1$$
then the continuous STFT is invertible. Several observations:

These are very mild conditions to satisfy. This should not be surprizing, the STFT is extremely redundant, so the information is everywhere in the time-frequency plane
yes, you need a second window in the synthesis time-frequency plane (but it can be the same)
The windows $h$ and $g$ are complementary in second-order. Many possible $g$ exists, and especially:
$g=h$, provided that $ \int_{-\infty}^{+\infty}h(t)h^*(t)\mathrm{d}t = 1$ which is easy if the window does not have zero energy.

In the discrete-time/discrete-frequency domain, this is far more complicated in the general case. Yet, for power-of-two frame length, for certain window (triangle, Hann) and overlap ($1/4$, $1/2$, $3/4$), there exist direct inverses, easily computed from the analysis window. One shall compensate a bit in amplitude (like in your case), by a factor related to overlap and redundancy. Those cases are often standard, as they are quite easy.
However, there are generalized inverses, that allow for less redundant decompositions, at a cheap prices.
Related answers:

Meaning of Hop-Size in Filter Bank Interpretation of Short Term Fourier Transform
Normalizing STFT output and calculating energy of signal


In MATLAB ,do the input signals need to be of the same sampling rate for the conv function to work correctly?

What an interesting question! An impulse response sequence, taken on its own, does not have a specific numerical sample rate (measured in Hz) associated with it. If I applied a [1, 0, 0, 0, ...] unit impulse sequence to a 3-point moving average lowpass filter then the filter's impulse response would be the sequence [1/3, 1/3, 1/3, 0, 0, 0, ...].  There's no specific numerical sample rate (measured in Hz) associated with either the input unit sample sequence or the filter's impulse response sequence. If you insist on talking about the somewhat unusual notion of the sample rate of a filter's impulse response I suppose you can say, "The sample rate of a filter's impulse response is always equal to the filter's input data sample rate."

For example, if to my 3-point lowpass moving average filter I apply an input sequence whose sample rate is 1000 samples/second (1 kHz), then the filter's output sequence has a sample rate of 1000 samples/second. On the other hand, if I apply an input sequence whose sample rate is 24000 samples/second (24 kHz) to the same lowpass filter, then the filter's output sequence will have a sample rate of 24000 samples/second. (WARNING: I suggest you never discuss the sample rate of a filter's impulse response during a job interview.)


The serial/parallel converter shown in most OFDM transmitter block diagrams is required in a hardware implementation but not very meaningful for software implementation. Because in software there is no IFFT block that actually works with parallel inputs, it's all serial processing in reality.

Let $S(k)$ be your input stream of complex numbers and $N$ be the number of subcarriers, i.e. the IFFT size. Then you have to calculate the IFFT of every $N$ samples in $S(k)$, i.e. of the vectors $[S(iN),\ldots,S((i+1)N - 1)]$, with symbol index $i=0,1,2,\ldots$. The result is the OFDM modulated time domain signal. If I understand your example correctly you have to apply the IFFT to both $S_1$ and $S_2$, yes.

The output of the IFFT operation is also complex, in general. Unless you use two audio channels followed by an Inphase/quadrature mixer that shifts the signal to a carrier frequency the vectors defined above have to exhibit a complex conjugate symmetry in order to obtain a real-valued time domain signal. (So-called discrete multitone (DMT))

The allocation of data to subcarriers happens through the IFFT algorithm you're using. The FFTW library, for example, allocates $S(0)$ to frequency zero and $S(N-1)$ to "frequency" $-N/2$. What physical frequency this corresponds to in turn depends on the sampling frequency of the digital-to-analog converter you're using. Orthogonality is inherent to the (I)FFT algorithm, you don't have to apply any further processing to obtain it.


Multiple different length Windows for high resolution realtime FFT

You can, and you don't need to zero pad, but a little overlapping does come into play.

Skipping the gnarly details and hitting this with a conceptual approach.  Applies to two or more "levels".

Start with a longer frame for your low tones.  The key is you need to smooth your signal fairly heavily.  Lots of ways to do this.  Then do a DFT (what a FFT really is) on the long interval.  The reason for the smoothing is to eliminate alias frequencies at this level.  The effects of your smoothing will be frequency dependent, so there are formulas for the different smoothing techniques.  Do the inverse of the FFT of the smoothed signal with attenuation compensation applied.  Subtract this from your original signal.  You have now documented and removed your lower frequencies.  Divide the interval into smaller chunks, and do DFTs as normal at the top level.

"Pasting" them together can also be done in different ways.  Easiest is just to display one spectogram above the other as "display bands".  You can also rescale the output from being linear in frequency to exponential if you want to put your results on an octave scale.


Hamming window in frequency domain and back

Your initial plot doesn’t show the full input to the FFT, which is 256 in length, and IIRC zero-padded.


How $\Delta f=\frac{1}{2T}$ satifies the orthogonality condition?

You're absolutely right, it must be $\Delta f = 1/T$, since this is the base frequency of the fourier transform of periodic functions with period $T$.

Otherwise, it's clear for $m=3$, for example, that the exponential would become

$g(t)=\exp\left(j \frac{2\pi}{2T}\cdot 3 \cdot t\right)$

So it does not obey periodic boundary conditions at the right boundary, i.e. at $t=k\cdot T$ because

$g(T)=\exp(3\pi j)=-1 \not= 1 = g(0)$

So the problem statement in the book is wrong, but seems to come to the right result with either wrong calculus or by already knowing the right result. Maybe the author confused $\Delta f$ with the Nyquist frequency or something. Or it was just a typo in the book.

Remark: the sinc function has zeroes for all integers $n,m$ except where $n=m$, which is kind of a pitfall of $sinc(x)$, because the limit for $x\to 0$ is one, not zero as the occurrence of the sine in the numerator alone would imply. So the orthogonality condition that you have mentioned in the question could be extended to an orthonormality condition

$$<s_{ml}(t),s_{nl}(t)>=\delta_{mn}$$

which actually also holds for $n=m$, where it presents the normalizability of the basis functions (requiring appropriate normalization factors, of course).


Why does an 8 point DFT behave differently from a 16 point?

Your sin frequency is 50 Hz, but in your 16 point DFT there is no 50 Hz bin.  There is a 0 Hz bin and a 100 Hz bin, and you can see that some of the signal power goes into both of them.  This is a non-ideal way of capturing a 50 Hz tone.

The 8 point DFT has a 50 Hz bin, and you can see that it captures almost all of the signal power.


Your understanding is correct and you can apply Hann windowing on discrete signals. The choice of your window depends on how easily you can visually resolve two sinusoids that are closely spaced in the spectrogram. $$STFT = X(m,\omega) = \sum_{n=-\infty}^{n=+\infty}x_nw_{n-mL}e^{-j\omega n}, $$ where $L$ is the step size of spectrogram (interval after which you place next window). and $W$ is the size of window, which decides resolution of frequency (row-wise separation of frequencies). Higher the $W$, more finer the resolution. You can see there is overlap between adjacent windows, the overlap being $W-L$.


Take a look at pixel art scaling algorithms, particularly the hqx algorithms.  They're typically limited to a integer scaling factor (i.e. 2x, 3x, 4x) so they aren't completely general purpose but they should do quite well on line drawings.


Is there a convolution mistake in my method?

You just multiplied the two functions and integrated them but you didn't convolve them. You must compute

$$y(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau=2\int_0^2h(t-\tau)d\tau-\int_2^3h(t-\tau)d\tau$$

Alternatively, you can compute the step response

$$a(t)=\int_{-\infty}^th(\tau)d\tau$$

Because $x(t)=2u(t)-3u(t-2)+u(t-3)$ the output can be written in terms of the step response:

$$y(t)=2a(t)-3a(t-2)+a(t-3)$$


How Can Telecommunication System Know Which Part Is Signal and Which Is Noise?

To elaborate on MM's answer just a little bit.  Unless the receiver has some sort of expectation of the nature of the signal there is no way to tell.  With an expectation (for instance a pure tone), only a best estimation can be made with the rest assumed to be noise.  With something like serial communication a voltage past a threshold is considered a "1" and no (or low) voltage is considered a "0".  In noisy channels, parity bits (or error correcting codes) are used to ensure "what was sent".  There are many possibilities, but the receiver has to know the nature of is being sent, which seems to be the answer you are fishing for.


I recommend a standard dq-PLL, check this reference https://vbn.aau.dk/ws/portalfiles/portal/273236528/PLL_Review_RSER.pdf

There are 2 required integrators. For algebraic loop reasons, the last integrator, the one that converts the angular frequency to the phase, should be a forward-Euler integrator. The integrator in the PI controller, can be a trapezoidal integrator if you want.
Edit : Ts = Sampling Period


How Discrete sequence could have more than one Frequency

When you are working with discrete sequences, frequency really becomes a different concept that "real world" frequency.  You can use analogy to relate information from discrete analysis back to the real world, but the two worlds are different.

First, once you have converted a real time / continuous signal to a discrete version, you immediately change the data representation such that no frequencies above 1/2 the sampling rate can be represented.

Second, the concept of frequency in the discrete world, with regard to the FFT, is based on the Discrete Fourier Transform (DFT).  The DFT imposes a restriction on your analysis.  The DFT assumes that all data is periodic.  The longest period represented in the transform is based on the length of samples being analyzed.  The representation says that your data is composed of a finite number of evenly spaced (in frequency) periodic sequences.  The set of periodic sequences are thought of as frequency bins.  The number of frequency bins is equal to to the length of the data.  So if you have say 8 samples your are working with, there are 8 frequency bins.  Half of the frequency bins are "reflections" of each other (represent the same frequency) so you are now down to 4 bins.  Two bins fall into the DC component, so you are left with 3 periodic sequences to represent your data.  So in terms of the original random signal you are analyzing, by doing an FFT analysis, you converted the signal from random to a band limited periodic signal composed of a finite number of frequencies.  This is very different from a real world random signal.

In your specific case, where you have computed fft([40 2 8 31 6 8 -4]), I don't know the mechanics of the specific FFT algorithm, but it looks like the data has been padded up to 8 samples, so you have a DC component and 3 additional periodic sequences.  Notice that the values ranging from 2 to 7 are symmetric around 4.5, so your periodic sequences are for frequencies represented by bins 2, 3 and 4 (5 is a reflection of 4, 6 is a reflection of 3 and 7 is a reflection of 2).  Bin 1 would be considered the period of the overall sequence.

You have to consider the sampling rate to convert the bins back to real frequencies:

bin# * SampleRate / NumberOfBins


How do we impose Hermitian Symmetry in optical OFDM?

I don't understand your question about the IFFT step: as in "normal" OFDM, you send a full vector to the IFFT, as that is a vector-valued operation of vectors. The output are (interpreted as) time-domain samples.
In fact, you really shouldn't have put Step 3 and 4 as separate operations, and certainly not in that order: the input to the IFFT has to be hermitian symmetric to yield a real-valued signal (which is what you need if you don't want to implement a quadrature upconverter to mix your signal into a carrier frequency, so in all honesty, fig. 3 is misleading in that the $-\sin(2\pi f_c t)$ branch doesn't have anything to do if you ensure your IFFT input is hermitian symmetric.)

My confusion here is if we have to "impose" that certain carriers must take a given value (e.g. carrier 0 and carrier N/2 must equal 0) then isn't that somehow "overwriting" the value for that carrier that we have applied at the QAM modulation stage?

No, you're just mixing up what you're doing: you're not free to assign an arbitrary value to each subcarrier at will, you need to make sure the vector is hermitian symmetric. You hence only have half the subcarriers to choose freely, the other one is inherently set to the complex conjugate.
So, you're not "overwriting" anything. You're just only assigning your QAM symbols to half the subcarriers.


Yes, what you've observed is correct:

The $N$-point DFT only "compares" a signal to the set of complex oscillations $e^{j2\pi f n}$ for frequencies $f$ that fit an integer number of times into the DFT length $N$.

Any other frequency will not be mapped "sharply" to a single DFT bin, but will have some energy in the surrounding bins and spectral repetitions.

We use windows on the input signal of a DFT to shape that "bleeding" a bit.

It is very common to increase the "resolution" by having a greater $N$. To do that, you simply take $K<N$ samples of your input, and attach $N-K$ zeros, so that you end up with a vector of $N$ samples. That way, you can have a finer "grid" of frequencies.


Getting Frequencies Corresponding to Peaks In FFT Plot - MATLAB

MATLAB's indexing begins at 1, instead of 0, as most programming languages use.  This causes a little wrinkle in interpreting the fft results.  Index 1 corresponds to bin 0, index 2 to bin 1, etc.  So to get the bin number, subtract 1 from the index.

The bin number corresponds to the frequency in units of cycles per frame, where frame is the section of the signal your are analysing.  Usually, the variable $N$ is used to denote the number of samples in the frame.  You also need to know the sampling rate, usually denoted by $f_s$ which usually has units of samples per second.  CD quality audio has a sampling frequency of 44,100 samples per second.

To calculate the frequency in hertz (Hz) of a bin number, usually denoted with a $k$, do the following calculation:

$$ f_{HZ} = k \cdot \frac{f_s}{N} $$

The units for the equation goes like this:

$$ \frac{cycles}{second} =  \frac{cycles}{frame}   \cdot  \frac{ \frac{samples}{second} }{ \frac{samples}{frame} } $$

Hope this helps,

Ced


You do not need any additional bits to help the decoder converge.  There should be no need for convergence at the boundaries, because the start and end states are completely known (the all zeros state) so there is no guessing involved.

If you are getting more errors at the end it is likely due to a bug in the decoder.


How to choose one among different types of LPF?

There are 4 main types of LP filters, Butterworth, Chebyshev1, Chebyshev2, and Elliptic filter. Each differs from the other in the oscillations it has, and in which part of the spectre it is.

Butter has the worst amplitude cutoff, eg the widest transition region, but no osclillations;
Cheb1 has better transition range, but has oscillations in the pass region;
Cheb2 is the opposite of Cheb1, has oscillations in the reject region, but not in the pass region;
Elliptic has the narrowest transition, so therefore the best, but has oscillations in both the pass and no-pass zone.

Heres a picture to help you visualize:

http://i.stack.imgur.com/rJMQw.png

Now, seeing as you need a peek detector, you should use a flat filter for the pass zone, so Butter or Cheb2. I advise Cheb2, because you don't really care about the higher frequencies.

In matlab, the function is cheb2ap, and is used like this:

[z, p, k] = cheb2ap(N, La)
% where La is the allowed ripple in the stopband (usually around 0-1.5 dB).

% And now the filter transfer function is acquired with:
b = k*poly(z);
a = poly(p); % b is the polynomial for zeroes, a is polynomial for poles
[H w] = freqz(b,a, N)  % N is the number of points for calculation

Now you have frequency response of your filter, H, and vector for its frequency, w, and you can plot the filter with plot(w, abs(H)).

Hope this helps :)


ROC of the function in the problem 9.14 of Oppenheim's Signals and Systems textbook

You're right, the ROC of the Laplace transform of a two-sided signal is a strip in the complex plane. In your case, the imaginary axis is inside the ROC, and the ROC is limited by the poles in the right and left half-planes. If the ROC were the right half-plane, the signal would be right-sided, which is clearly not the case.


How to get rid of discontinuities at the start of a signal caused by the filter

If you are post processing you can instead use a "zero-phase" filter which passes the signal through the filter twice in the forward and reverse direction, resulting in an output that is completely aligned with the input as would be expected with a filter with no delay, and with a significant reduction in start-up deviations (this is a non-causal filter hence the restriction on postprocessing).
Such a filter is available in MATLAB, Octave and Python scipy.signal as filtfilt which has further documentation on its use.
I demonstrate this as shown with the plot below using a traditional 91 tap FIR low pass filter for an in-band cosine signal and again with the same coefficients with filfilt (in red). Notice there are no transients or delay.

As far as any residual transient as @MBaz astutely points out, below is the plot zoomed in on the first 4 samples showing the filtered output as (+) and the original signal as (o). The line plot simply connects the output samples. This is due to the fractional cycle remaining at the end of the sampled data where I did worst case with a 1/2 cycle residual in the capture; if done with an integer number of cycles there would be no such mismatch.

And the same portion of the waveform later in the sequence showing close alignment.

However in comparison to the blue line in the first plot even in this worst case condition, we do see the significant reduction in the start-up transient due to the use of the zero-phase filter (limited to post-processing applications). The alignment of the input and output is also a huge plus for analysis applications where we many otherwise have to go through additional efforts to offset our waveforms for the delay when comparing input to output.


wavwrite with 32 bits per sample seems to overflow?

It seems like it is a bug. Octave's wavwrite() function writes 32-bit wav in 32 bits int (type 1 format) instead of normalized 32-bit floats (type 3 format). Unfortunately, 1 can't be represented in 32-bit type 1 format which is why you get this result.
The solution could be to use the audiowrite() function, or to write 0x7fffffff instead of 1.
More about wav formats.
EDIT: seems from the bug report that the audiowrite function has the same problem...


Real and Imaginary data from Network Analyzer

The S paramater data are complex numbers that have magnitude and phase, or stated another way, real and imaginary components. The way we typically see the data presented on a network analyzer is with a Magnitude vs Frequency plot, and a Phase vs Frequency plot. If we capture the data we may see it listed as Real and Imaginary. The two are related as follows:

Given: Magnitude (K) and Phase ($\phi$)

$Sxx(\omega) = K(\omega)e^{j\phi(\omega)}$

(showing magnitude K vs frequency together with phase $\omega$ vs frequency since K at angle $\phi$ is $Ke^{j\phi}$ in compact exponential form.)

Has real and imaginary components that are also functions of frequency as follows:

$K(\omega)e^{j\phi(\omega)} = I(\omega) + jQ(\omega)$

Where I is often used for the real component (In- phase), and Q is often used for the imaginary component (Quadrature phase).

The relationship between the two is given by Euler's identity:

$K(\omega)e^{j\phi(\omega)} = K(\omega)cos(\theta(\omega)) + jK(\omega)sin(\theta(\omega))$

The real and imaginary numbers define the magnitude and phase versus frequency for the relevant S parameters. For example S11 defines the reflected signal with values between 0 and 1 (for passive devices) and phases anywhere from 0 to 360 degrees. An open circuit for example will reflect all incident power back toward the source in phase (so would be 1 (0 dB) at angle 0 degrees), while a short will reflect all the incident power back toward the source out of phase, 180 degrees (so would also be 1 but at angle 180 degrees).  S21 is the transmission from port 1 to port 2 of a 2 port network, and a simple cable is a good demonstration of phase: a lossless cable is just a delay, and a delay has a linear phase with frequency so would be magnitude 1 (0 dB) with a phase that increases linearly negative as frequency increases.

Thus to demonstrate the significance of the phase, you could predict the magnitude response of a comb filter that is created with two power dividers and a long length of cable between one leg and short on the other (split, go through the cables and then combine for a resulting two post network), by measuring S21 phase through each of the cables.

At certain frequencies the phase between the cables will be 180 degrees and thus S21 magnitude through the entire network will be nulled when the two signals combine, while other frequencies they will be in phase and the signal will be combined, thus resulting in a comb pattern in magnitude such as that shown in the normalized frequency plot below.

The phase of the reflected signal (S11) and transmitted signal (S21) is used often in many aspects of RF design including impedance matching and filtering design (since like in digital filter design filters can be constructed through the summation of weighted delays, and hence weighted phase shifts). Having a deep understanding of this is fundamental to RF design.


Trying to decompose a signal into sine waves whose wavelengths aren't restricted to 1/N times sample range?

There is absolutely no need for explicit interpolation. The DFT can do that for you. Take you signal frame and apply a windowing function to minimise edge artefacts. Afterwards pad the frame with zeros to determine the number to output bins of the DFT. Then apply your FFT and you have an interpolated spectrum with the frequency step size you desire.


Why is fourier analysis used in doppler effect?

The fourier transform of a signal (like sound) will show you the sound in the frequency domain. In other words, the independent variable, or x axis, will be in terms of frequency (such as hertz or radians per second), and the y axis will show the magnitude of those frequencies.
One reason that the fourier transform is useful is that it will show what frequencies are present in a signal and in what quantities. When used in the experiment you described (if I understand it correctly), you could help demonstrate the doppler effect by showing that the sound source is outputting a tone of a constant pitch, and in the recording of the sound source moving in the circular pattern, the fourier transform of the recording will show that the frequency is not constant, and appears in significant quantities in the sound.
If you record the sound as the source moves in a circle, then the fourier transform of the sound should appear to be high in a range close to the pitch that the source is outputting.
A tool I frequently use for computing the fourier transform of an audio file is matlab. I am sure there are many tools online to help you compute this as well if you do not have access to matlab.


Your vector

t = linspace(0, T, num_samples);

is the problem. It should be

t = (0:num_samples-1)/fs;

Note that, even though you claim it in a comment in the code, in the current version, dt = t(2) - t(1); is not equal to 1/f_s.


Equal contribution criteria during image pyramid construction

The following picture from this set of lecture slides explains the equal contribution. The contribution from two nodes in the same layer must be equal. With reference to the figure, the sum of weights on green edges must match the sum of weights on the blue edges.


Yes it does exist and these are called "all-pass" filters in that over the band of interest the magnitude does not change (other than a possible fixed gain at all frequencies) but modify the phase response. The time delay at a specific frequency is the negative derivative of phase with respect to frequency, so the goal of these filters is to selectively change the slope of the phase according to the desired delay at each frequency.

The primary characteristic of a digital all-pass filter is that it will only have poles inside the unit circle on the z-plane and for each of those poles will have a conjugate reciprocal zero outside the unit circle, as depicted in the plot below.


Can Oversampling help to mitigate the effect of additive white gaussian noise?

Oversampling will not do anything to change the existing noise in the waveform other than the quantization noise that is added. Sampling does not change the waveform but the quantization process of sampling with ADC’s will add additional noise to it (quantization noise and other additive noise collectively given as the “Equivalent Number of Bits” or ENOB). The OP is referring to the Gaussian noise that may be sampled, which itself will not be modified other than the additional noise added. This makes complete sense since our goal with sampling a waveform is to represent it with minimum distortion; there is nothing different between noise and a waveform of interest that would somehow be sampled differently by the ADC.
With regards to the quantization noise, oversampling will spread the additive quantization noise out over a wider frequency and thus reducing it density on a per Hz basis, since the quantization noise is well modelled as a white noise over the sampling bandwidth. The total quantization noise power doesn't change, but the noise density would be lower since that noise is spread over a wider bandwidth. When the ENOB is limited by spurs, then further oversampling will not reduce the noise further since a narrow band spur is confined in frequency rather than a spread noise. Further as sampling is increased the ENOB typically reduces due to non-linearities and aperture uncertainty in the ADC front end.
Given this, the typical goal in a system is to have sufficient gain in front of the ADC such that the amplified thermal noise floor is higher than the added noise from the ADC with sufficient margin such that the noise floor we observe digitally is the thermal noise floor of the analog waveform and not the added quantization noise of the ADC (for the latter case we would then be swamping out the signal with additional added noise which is typically not desirable when our interest is SNR). Oversampling can help meet this goal by effectively increasing the number of bits in the ADC (increasing its dynamic range).
For further details please refer to: What are advantages of having higher sampling rate of a signal?


Bilateral $\mathcal Z$-transform of exponential

In complete analogy with the bilateral Laplace transform of $x(t)=e^{-at}$ (which doesn't exist), the bilateral $\mathcal{Z}$-transform of $a^n$ doesn't exist either. The series

$$\sum_{n=-\infty}^{\infty}a^nz^{-n}$$

converges nowhere, simply because $a^n$ grows without bounds for $n\rightarrow -\infty$ if $|a|<1$, or for $n\rightarrow\infty$ if $|a|>1$. Of course, for $|a|=1$ there series doesn't converge either.

EDIT:

As for your computation of the $\mathcal{Z}$-transform of $a^n$, the mistake lies in the fact that in addition to the algebraic expression of the transform you also need to consider the region of convergence. If you split $a^n$ (as you did) as

$$a^n=a^nu[n]+a^nu[-n-1]\tag{1}$$

you can compute the $\mathcal{Z}$-transform of both right-hand side expressions separately:

$$\mathcal{Z}\{a^nu[n]\}=\frac{z}{z-a},\quad |z|>|a|\\
\mathcal{Z}\{a^nu[-n-1]\}=-\frac{z}{z-a},\quad |z|<|a|\tag{2}$$

Note that the region of convergence (ROC) for the first part is outside the circle with radius $|a|$, whereas the ROC of the second part is inside the circle with radius $|a|$. The ROC of the total expression would be the overlap of the two ROCs, which is zero. Consequently, the sum doesn't converge anywhere and the $\mathcal{Z}$-transform of the total expression doesn't exist.


Implementing Wavelet Transform using Equations

I found the Wavelets for Computer Graphics: A Primer a good introduction to the Haar wavelet and its role in image processing.


How to expand on and deepen a basic skill set for image/video processing?

I hope you realize (from the form of your question, I think you do) that three months of reading up about it doesn't really make you an image processing expert.

I am familiar with topcoder, but while the difficulty of the problems and approaches that you needed to use there might be similar to those of image processing problems, to design a usable image processing application: do problem specific research, implement your findings and new ideas, test it to get reliable results, you need much more than a couple of hours (or even a week as for the topcoder marathon match).

If you want to just learn more about the basic image processing tools I suggest you:

take all your literature, random-open a page and study and implement the concept explained there
take a computer vision library such as OpenCV, or work in matlab, take some images, and try using and displaying various image processing methods on them (from simple morphology to more complicated things like Hough line detection). Try and predict results, graphically display them, and understand what you get.

If you are interested in specific problem in image processing (e.g. I am doing Content based image retrieval - CBIR right now -- and by that I mean, almost for the last 5 months), and want to get better in this specific topic, then my advice would be something like:

search google scholar, Mendeley article database, IEEEXplore for relevant works on the subject
download way too many papers, maybe 3-5 as many as you are prepared to read. Skim abstract, intruduction and conclusion sections and filter out the articles that sound relevant and interesting. Give priority to current (recent) articles, as well as articles with current references.
read all of that. Read the most important ones more than once. The one that will be the basis of your work is expected to have margins filled with your scribblings and to look at least a little bit scrappy by the time you finish your 5th read-through.
implement that. Go from simple, straightforward implementation to the optimized implementation.
before you start doing this, ensure that you have a dataset for testing. Testing is not relevant if it is not done on a large enough dataset. Some works mention the datasets they are using, or their datasets can be found on affiliated university websites.
there are public datasets (e.g. this one for object classification) that also organize annual challenges where you can test your new ideas against many peoples new ideas (but that's way bigger than topcoder :D)
if you get decent results, that's good. If you get better results than state-of-the art using some of your new ideas, double check them. Then triple check them. And then publish a fancy article ;)

I'm sure you can find a middle between just learning something about image processing just for fun and doing real time-consuming research from what I wrote... Actually, here's an idea: hang around here and try to understand and help to solve other people's problems! All of them need research-think-apply steps, they just vary with the depth of the steps they require ;) In any case, I hope it helps.


how can i perform a registration to 2D images using ''rigid " transformation in matlab?

A rigid, or isometrie transformation, is a transformation in which the distances between the points that are transformed, are preserved. This means that an object formed by these points will have the same shape and size after transformation. Examples are rotations and translations.

A 2D column vector can be rotated by $\theta$ as follows:

$\left(\begin{array}{c}x' \\ y'\end{array}\right) = \left ( \begin{array}{cc} \cos \theta & - \sin \theta \\ \sin \theta & \cos \theta \end{array} \right ) \left(\begin{array}{c}x \\ y\end{array}\right)$

In your case you will probably be seeking a combination of rotation and translation between two images meaning that the relation between the points of these images can be defined as follows:

$\left(\begin{array}{c}x' \\ y' \\ 1\end{array}\right) = \left ( \begin{array}{ccc} \cos \theta & - \sin \theta & t_x\\ \sin \theta & \cos \theta & t_y \\ 0 & 0 & 1 \end{array} \right ) \left(\begin{array}{c}x \\ y \\1\end{array}\right)$

Or short:

$\textbf{x}' = \textbf{A}\textbf{x}$

This particular transformation has 3 degrees of freedom, one for rotation and two for translation. This means that you can solve this matrix for 2 point correspondences $\{ \textbf{x}' \leftrightarrow \textbf{x} \}$. That means in your case that you need two points in the original image and two matching points in the distorted image to ascertain the translation and rotation between them.

Ultimate reference: Multiple View Geometry

EDIT

The deformation between these images is obviously not an isometrie. For one the branches appear and disappear between them. I played a little with a slightly more involved form of feature based registration based on the assumption that the distortion is a similarity transform with 4 degrees of freedom, and the features will make up for the difference in edges. In short:

SIFT feature detection and matching
RANSAC estimation to detect inliers
Direct Linear Transform (DLT) variant for estimation of scale and rotation
Run of the mill cross correlation to detect the translation

Results, not brilliant, but here goes

All code on Github, demo file is so_13_02_2012.m.


When using a moving average filter, or rectangular convolution in the time domain, the time domain smoothing and so-called spectral leakage are nearly opposing effects.  The smoothing is due to removing high frequency content, and the leakage is the name given for portions of this frequency content that end up not being removed.  Since sharp peaks and transients can be partially composed of high frequency content, smoothing by this kind of convolution can diminish them.

Time-domain convolution windows other than rectangular may remove less of the very nearby frequencies, but better remove more of the far away frequency content (relative to some transition frequency, which is related approximately to the window width).

An adaptive polynomial regression will remove different spectral characteristics from a signal at different points in time, which could either expose or hide "important" elements of the signal, depending on the signal, and what you are looking for.


How does averaging increase the Signal to Noise ratio?

Intuitively this is true, because averaging a zero mean noise processes approximates its expectation value - which is zero. More rigorously:

If the signal $x$ that you want to observe (estimate, actually) is constant for all observations $y$ we can write the $n\text{th}$ observation as
$$
y_n = x + r_n
$$
where $r_n$ is the noise which is different for every observation, of course. The SNR $\gamma_1$ of $y_n$ is
$$
\gamma_1=\frac{\operatorname{E}[|x|^2]}{\operatorname{E}[|r_n|^2]}
$$
where $\operatorname{E}$ is the expectation operator. By averaging $N$ observations of $y_n$ we get
$$
\hat x = \frac{1}{N}\sum_{n=1}^N( x+r_n) = x + \frac{1}{N}\sum_{n=1}^N r_n
$$
The SNR $\gamma_2$ of the averaged observations is
$$
\gamma_2=\frac{\operatorname{E}[|x|^2]}{\frac{1}{N^2}\operatorname{E}[|\sum_{n=1}^N r_n|^2]}
$$
To simplify that expression, we assume that the noise sequences $r_n$ are uncorrelated (1). In that case, the expectation of the squared sum of $r_n$ is equal to the sum of the expectations of the squared $r_n$, i.e.
$$
\operatorname{E}\left[|\sum_{n=1}^N r_n|^2\right] = \sum_{n=1}^N \operatorname{E}\left[|r_n|^2\right]=N\operatorname{E}\left[|r_n|^2\right]
$$
Inserting this in the expression for $\gamma_2$ yields
$$
\gamma_2=N \frac{\operatorname{E}[|x|^2]}{\operatorname{E}[|r_n|^2]}
$$
Comparing the SNR $\gamma_1$ of the single observation and the SNR $\gamma_2$ of the averaged observation we see that averaging increases the SNR by a factor of $N$.

(1) This holds for a Gaussian noise process, for example. Right now I'm not sure for what classes of stochastic processes it holds. Comments on this are welcome.