Partial Fraction Expansion for Inverse Fourier Transform

You should indeed replace $e^{j\omega}$ by the more general complex variable $z$ because when talking about poles and zeros of a rational function, you must consider the whole complex plane. The fact that $e^{-j\omega}=3$ doesn't have a real solution just means that the corresponding pole does not lie on the unit circle.

After having performed the partial fraction expansion, there is nothing that keeps you from setting $z=e^{j\omega}$.


Hopefully you understand why and in what situations those methods are used.

The integer 'K' in the Bartlett and Welch methods is not arbitrary.  The value of 'K' is a "trade off" that you must make. In both methods, the greater is 'K" the lower is the variance in your final DFT results. That's a good thing. What the PSD's author doesn't tell you, though, is that the greater is 'K" the larger is your final DFT result's frequency-domain sample-to-sample spacing (measured in Hz). And that's usually a bad thing. So you must decide what value of 'K' is acceptable to you. I'm not qualified to comment on the Blackman-Tukey method.


Linear Predictive Coding - Meaning of spectrally Flat

1) what does spectrally flat mean ?

Spectrally flat means that the signal has energy at most / all frequencies, and this energy is roughly the same level across the frequency range.

  2) why does an impulse train have a flat spectral ?

As you're probably aware, the Fourier transform of an impulse train is also an impulse train. That means that the claim by the link that the impulse train is spectrally flat is not correct.

The reason an impulse train is used in modeling speech is because it gives a tonal quality for non-fricative sounds (mostly vowel sounds). That tonal quality is hard to generate by filtering if the source being filtered is white noise.


CZT: Do Bluestein algorithm perform backward transformation?

It is easy to calculate IDFT using DFT, by applying the property:
$$\operatorname{IDFT}(x) = \frac{1}{N}\big(\operatorname{DFT}(x^*)\big)^*,\tag{1}$$
where $^*$ denotes elementwise complex conjugation. You can find this equation in a more verbose form in Eq. 13-55 of Richard G. Lyons "Understanding Digital Signal Processing" 3rd edition.
Probably this was used along with conjugation properties of the exponential.


Pole Magnitude and Damping Ratio relationship

The magnitude of the poles alone cannot determine the damping ratio, which in turn is related to the rate at which the oscillations about the steady state value will decay.
In fact it is the ratio of the absolute value of the real part of the poles and the magnitude of the poles that gives the value of $\zeta$.
If you think about a second order system its characteristic polynomial would be:
$$
s^2+2\,\zeta\,\omega_n\cdot s \,+\, \omega_n^2 = \\\left(s\,+\omega_n\cdot \zeta-j\,\omega_n\,\sqrt{1-\zeta^2}\right)\cdot \left(s\,+\omega_n\cdot \zeta+j\,\omega_n\,\sqrt{1-\zeta^2}\right)
$$
with poles: $p_{1,2}\,=\,-\omega_n\cdot \zeta \pm j\,\omega_n\,\sqrt{1-\zeta^2}$
However, the above pole forms are valid for $0 \le \zeta \lt 1$. These systems are called underdamped.
Whereas the damping ratio is:
$$
\zeta = \dfrac{|\Re(p)|}{\omega_n}
$$
As a result of the above the magnitude and the damping ratio relate to each other through:
$$
\zeta = \dfrac{|\Re(p)|}{||p||}
$$
You are correct in that the angle of the pole, with respect to the left hand side real axis relates to $\zeta$, and this is the formula:
$$
\cos{\theta} = \dfrac{|\Re{p}|}{\omega_n} = \dfrac{\zeta\cdot\omega_n}{\omega_n} = \zeta
$$
The conclusion that the magnitude alone, whether small or big, cannot determine the decay rate of the oscillations is also apparent from the fact that if we normalize the poles with the natural frequency $\omega_n$ their magnitude will be always equal to 1, this means that if the angle of the pole somehow varies, then they will always be on the unit circle. In such case as the angle varies the damping ratio will also vary but the magnitude will always be the same no matter what the angle.
So its the real part of the normalized pole that determines the damping ratio
rather than its magnitude, alone.
Note: A slow, from a steady state point of view, pole is in fact one that has a small real part, while a fast is one that has larger real part, in absolute value.
Systems with $\zeta = 1$ are critically damped, having conjugate imaginary poles with magnitudes at the system's natural frequency:
$$
p_{1,2} = \pm\, j\, \omega_n
$$
So this means that the hypothetical system has a pair of conjugate poles with equal magnitude:
$$
||p_{1,2}|| = \omega_n
$$
Those systems may exhibit oscillatory behavior sustained longer than a underdamped system, for example, a sinusoidal response. The poles angle is $90^o$ degrees because their real part is zero. Hence $\zeta = \cos(90^o) = 0$ implying an undamped system but the magnitude of the poles is $\omega_n$, thus it depends on the natural frequency of the system which could be huge in value.
In many engineering applications, the measure that is used extensively for pole quality characterization is the $Q$-factor.
In fact, the quality factor and the damping ratio relate to each other through the following formula
$$
Q = \dfrac{1}{2\zeta}
$$
In addition, you may come across transfer functions with characteristic polynomials parameterized with $Q$ instead of $\zeta$.
$$
s^2+2\,\zeta\,\omega_n\cdot s \,+\, \omega_n^2 =
s^2+\dfrac{\omega_n}{Q}\cdot s \,+\, \omega_n^2
$$


Does any error correction code still work in such situation?

I'm afraid that this will be an incomplete answer, but I'm giving it in the hopes that it is better than nothing.

It is impossible to recover any information with a BER of 50% because, as information theory tells us, no information is getting through such a channel.  You can see this intuitively if you realize that no matter what bit you send, 0 or 1, it has equal chance of arriving at the receiver as a 0 or 1, so what you receive tells you nothing about what was actually sent.  In such a case it is theoretically impossible to create an ECC that will help at all.

Erasures are, in some situations, better than errors though.  Reed-Solomon decoders often use erasures to increase their error correction capability.  The reason why it works is because with RS codes you have two sets of unknowns: the error locations and the error magnitudes (the number of errors is also an implicit unknown).  Erasures give the decoder the error location so it only needs to figure out the magnitude.

Talking about "erasures" in the context of bits, though, makes little sense because with bit-based codes there is only one set of unknowns: the error locations.  The error magnitudes are always 1.  So if you have a bit "erasure" you are trading one unknown (error location) for another (error magnitude).

With bit-based codes like Viterbi codes the analog to erasures is using soft-decisions.  Soft-decisions contain "confidence" information- i.e. how sure the receiver is that the bit is really a 1 or really a 0.  That information can be used to greatly improve decoder performance.

All that being said, though, if your error/erasure rate is really that high you will probably not get good performance out of any ECC.


How to visualize this statement regarding Conjugate Symmetry

For $(1)$, how can I visualize this to see why the negative frequency components do not matter anymore$?$, does this mean the Fourier series representation of $x(t)$ becomes :
$$
x(t)=\sum_{k=0}^{\infty}c_{k}e^{jk\omega_{0} t}\;???
$$

No, the Fourier series is still what it was before this conjugate symmetry nonsense reared its ugly head. But what happens is that
\begin{align}
x(t) &= \sum_{k=-\infty}^{\infty}c_{k}e^{jk\omega_{0} t} & {\scriptstyle\text{original definition}}\\
&= \sum_{k=-\infty}^{-1}c_{k}e^{jk\omega_{0} t} + c_0 + \sum_{k=1}^{\infty}c_{k}e^{jk\omega_{0} t}& {\scriptstyle\text{breaking up the sum}}\\
&= \sum_{k=1}^{\infty}c_{-k}e^{-jk\omega_{0} t} + c_0 + \sum_{k=1}^{\infty}c_{k}e^{jk\omega_{0} t}& {\scriptstyle\text{changing indices}}\\
&= c_0 + \sum_{k=1}^{\infty}c_{k}e^{jk\omega_{0} t} + c_{-k}e^{-jk\omega_{0} t} & {\scriptstyle\text{combining sums}}\\
&= c_0 + \sum_{k=1}^{\infty}c_{k}e^{jk\omega_{0} t} + c_{k}^*e^{-jk\omega_{0} t} & {\scriptstyle c_k = c_{-k}^* \implies c_{-k} = c_k^*}\\
&= c_0 + \sum_{k=1}^{\infty}c_{k}e^{jk\omega_{0} t} + \left(c_{k}e^{jk\omega_{0} t}\right)^* & {\scriptstyle a^*b^* = (ab)^*}\\
&= c_0 + 2\sum_{k=1}^{\infty}\operatorname{Re}\left(c_{k}e^{jk\omega_{0} t}\right)
\end{align}
showing that one needs to compute or know the $c_k$ values for nonnegative $k$ only; the values of $c_k$ for $k < 0$ are not needed. Anal-retentives who believe in crossing their eyes and dotting their teas can always determine the values of the $c_{-k}$ by taking the complex conjugates of the known $c_k$'s.


At what time stretching ratio would an audiophile notice a wrong pitch

Audiophiles don't hear anything unless they are told what to expect.

Joke aside, this is equivalent to transposition by a bit less than a semitone. The effect won't be noticed on voices - this is way too low for any kind "chipmunk" effect to be observed. As for music, some people very familiar with the original music might detect the change of tonality and/or tempo. But such stretching ratios are not unusual at all - DJs doing beat-matching, TV channels speeding up a sitcom to squeeze in more ads, and many, many movie to video conversions (telecine) have routinely used it.


Is there a systematic method for converting an even length FIR filter to odd length?

Standard even-length filters are somewhat inherent to 2-band real FIR orthogonal wavelets. And they are most often non symmetric.
If you really want odd length, or at least central symmetry, most probably you will have to relieve one of the above assumptions. The simplest could be to use some biorthogonal wavelets. The classical JPEG2000 9/7 pair is probably close enough to orthogonality for most purposes. Other filters can be tabulated, see for instance: Wavelet Biorthogonal 2.8 (bior2.8).
Then, there is a possibly to have complex filters, but honestly I am not familiar anymore with the literature.
You can also use $(2M+1)$-wavelet filter-banks, to quit the dyadic case. There are a couple of $3$-band wavelet designs with odd-length filters.


Condition for optimal value to be the average

You don't need to assume anything about the data in order for the average to be optimal (in a least squares sense) for certain estimation problems. As an example, assume you want to estimate a constant value $x$ and you are given $N$ noisy measurements:

$$y_n=x+e_n,\quad n=0,1,\ldots,N-1\tag{1}$$

where $e_n$ is the error of the $n^{th}$ measurement. If our goal is to minimize the mean squared error

$$\frac{1}{N}\sum_{n=0}^{N-1}e_n^2=\frac{1}{N}\sum_{n=0}^{N-1}(y_n-x)^2\tag{2}$$

then the optimal estimate $\hat{x}$ is obtained by minimizing (2). Taking the derivative of the right-hand side of (2) with respect to $x$ and setting it to zero gives

$$\frac{2}{N}\sum_{n=0}^{N-1}(y_n-\hat{x})=0\tag{3}$$

From (3) the optimum estimate is obtained as

$$\hat{x}=\frac{1}{N}\sum_{n=0}^{N-1}y_n\tag{4}$$

which is simply the average of all observations.


Let us consider two random variables $X$ and $Y$. You are aware of the following terminologies with respect to (wrt) the random variables:

(Zero Mean) : $\mathbb{E}(X) = 0$ and $\mathbb{E}(Y) = 0$,

(Uncorrelated) : $\mathbb{E}(X Y) = \mathbb{E}(X) \mathbb{E}(Y)$. It is interesting to figure out how it is different from independence.

Now, an image specified, say $g(x,y) = f(x,y) + n(x,y)$ is a discrete random random process. The notion of zero mean and uncorrelated are not very different here, as for any two points $(x,y) = (x_1,y_1)$ and $(x,y) = (x_2,y_2)$, the image intensity at these pixel locations are random variables.


Maximum likelihood estimator for multiplicative Gaussian noise

OK, let's have a look at one of the problematic terms:
$$
\frac{\delta}{\delta x} \bigg[ \bigg(\frac{\tilde{d}[n]}{x}-\mu_A\bigg)^2 \bigg ]
= - \frac{2 \tilde{d}[n] \bigg (\tilde{d}[n] - \mu_A x\bigg) }{x^3}
$$
which can be verified by Wolfram Alpha.

The full derivative of the summation term is then just this summed over $n$.


Multiple parallel IIR bandpass filters with different centre frequencies combined into one filter for optimisation... possible?

So, what I took away from my multirate lectures is:

You don't use IIRs in that case, because the recursive part isn't shareable among multiple filter implementations (being specific to the output of each filter), and hence, you end up with a lot more multipliers than you'd have if you went for a FIR system and removed redundancies. This is especially true for hardware/FPGA implementation, where placing a wire to use a result twice is cheap. You're doing software, so saving a few registers worth of filter taps really won't pay -- you're breaking the linear nature of memory access that way, and that will impose a performance hit that is much, much worse.

In fact, I'd challenge you to look at how you can understand your filter bank as an array of filterbanks with identical filters.

For example, assume that all your BPFs have equally spaced center frequencies $\frac{f_\text{sample}}{N}$. Also, they might have different bandwidths, but since you say they are related to each other, let's assume there's one "maximum bandwidth" that is a multiple of all smaller bandwidths.

Now, implementing a polyphase filterbank ("channelizer") with a BPF of that maximum bandwidth and the constant frequency spacing is very efficient. Basically, it's like taking a single very BPF at full rate, add one FFT, and you get the complexity of the system for up to $N$ channels (note: this is $N-1$ free channels!). Then, take the channels whose bandwidth you'd need to be reduced further. Now, assume these $m<N$ channels all need to be reduced in bandwidth by a factor of $4$, for example. So that's $m$ identical quarterband lowpass filters! Implementing identical filters in hardware is pretty efficient, usually, especially because after the first filterbank, sampling rate might already be reduced (so you can multiplex multipliers); now, you're seemingly doing audio on a PC:

Take these $m$ channels, and just implement individual filters. There's nothing much you can or need to optimize here; assuming this is on something that has a PC-like CPU, caching will keep your filter taps close to CPU registers, and and only the actual samples will need to go through memory buses. Also, not that things like $\frac{1}{M}$ low pass filters are relatively relaxed in transition width requirements, and hence, FIR implementations are short in taps, and thus fast.

As usual in signal processing, decimate whenever possible! The polyphase filterbank might already decimate, but also your $m$ channel low pass filters should do that, too, and in a clever manner (a polyphase lowpass decimator runs at the output, not the input sampling rate!). Not processing samples you don't need is probably reason Nr. 1 that modern DSP is so capable.


Why is the error between the desired signal and estimated signal in the case of LMS filter remaining constant even after n number of iterations

As mentioned in the comment, I modified the code given here and was able to adapt the LMS filter with error tapering to zero. The only assumption I made is that (since I am not an audio expert and do not know how the channel from speaker to the microphone would look like), I assumed a 10 tap channel with only first 3 non-zero values (multi-path reflection from walls). OP is free to use a channel of his/her choice. Using this I generated 'actual_fb_path' by convolving white-noise with assumed channel. If I get access to 'actual_fb_path' it would be great. I would request OP to modify this as per the system model to see if it adapts to the delay.

clc
clear all
close all
training_input = rand(1,100);
h_actual = [0.8 0.3 0.05 zeros(1,7)]; %channel assumption
actual_fb_path = conv(training_input, h_actual);
d = actual_fb_path;
mu = 0.1;
%ha = lms(40,mu);
%[y,e] = filter(ha, training_input, d);
lms = dsp.LMSFilter;
lms.StepSize = mu;
lms.Length = 10;
[y,e]=lms([training_input zeros(1,9)]', d');
%subplot(2,1,1)
N=length(training_input) + length(h_actual)-1;
figure()
plot(1:N,d,'r',1:N,y,'b',1:N,e,'g')
title('System')
legend('Desired','Output','Error')


Inverse fourier transform of Hermitian function, getting an imaginary part

A (near) sinusoid in one domain (zoom in off center, and you might see it) can indicate a shift away from symmetry in the other domain.  So your initial imaginary signal may be 1 or 1/2 samples off from exact circular symmetry around element 0 (maybe around element 1 in MatLab or Fortran weird array indexing)


Questions about Active Noise Cancellation

That's a pretty broad question. Let's start with a noise cancelling head set, that's about the easiest device.

A ANC headset has an internal microphone, that's placed as close to the ear canal as possible. It than runs a regular control loop (http://en.wikipedia.org/wiki/Control_theory): It calculates a speaker signal that tries to keep the signal at the microphone at zero. If you want to play music at the same time, you simply make the "desired signal" to be the music and not zero.

However, there are all sorts of implementation problems

Designing the control loop is difficult in terms of stability, cancelling performance, spectral fidelity, etc.
Manufacturing tolerances
Control loops need to be fast, so they work better at low frequencies. Most good headsets use a well tuned combination of active and passive cancellation. Passive works a lot better at high frequencies
While the headset cavity is small, it still has a fair amount of acoustic resonance and wave propagation, so you need to correct for the sound at the microphone vs. the sound at the ear drum
Many newer headsets use outside microphones as well, so it's a combination of feed-forward and feedback cancellation

Once the spaces get larger than a headset, things get a lot more complicated. Now you need to take into account the spatially distributed sound field, the environment, all kinds of propagation modes.


How do I know if the given signal is discrete?

Fortunately, one typically does not receive signals from strangers in real life. One knows because one is working on an application.

A discrete signal can be equivalent to a continuous signal so you need a set of assumptions tying the two together, like a Nyquist sampling criteria and then you need to test for equivalence under those assumptions.

A function that only exists at discrete points( like time points) and there is no meaningful notion of intermediate values between those points is discrete. Some financial time series are like that.  Arrival times of a Poisson process is another discrete function.  If intermediate values on a continuous range are meaningful, a continuous function is likely.


Cropping binary images based on pixel counts in python

One approach that is always available for binary images are morphological operations.

The critical piece of code is
# from skimage import morphology
# import numpy as np

# filter elements smaller than 1/6-th of the total width
selem = morphology.rectangle(1, img_binary.shape[1] // 6)
img_opened = morphology.opening(img_binary, selem=selem)

hist = np.sum(img_opened, axis=1)
non_zero = np.argwhere(hist > 0)
top_bound = int(np.round(np.min(non_zero) * (1 - 0.05))) # 5% tolerance
bottom_bound = int(np.round(np.max(non_zero) * (1 + 0.05))) # 5% tolerance

Here is the notebook I used to produce the image above:
https://gist.github.com/FirefoxMetzger/5d0d5e207bf8182f61374ac9ad80f10f


Analog AGC is mainly used to increase the signal level to match the
ADC sensitivity? Is Analog AGC also used to attenuate the signal?.

For "Analog AGC" we use analog control components (typically voltage variable amplifiers and voltage variable attenuators) to adjust the receiver gain, and we use an analog power detector to measure the total band power (within the bandwidth of the filtering up to the detector). This measured power is compared to a target which generates an error signal (often with an op-amp in an analog circuit) and the error is integrated in a loop filter (also often with the same op-amp) which then integrates to the control voltage needed to make the error go to zero (and the loop filter can be further modified for higher order loop control). In the end the band power AS MEASURED AT THE DETECTOR will be constant except for dynamic conditions during power transitions.  Depending on what is used in the receiver for gain control, the AGC can amplify or attenuate all the signals in band as allowed.

Digital AGC is the first block used after ADC. Is the main purpose of
digital AGC is to provide constant level to the signal processing
receiver chain( RRC, CIC, demod....)

Digital AGC is not necessarily the first block after the ADC, or at least the level detector part of the AGC loop. A significant benefit of doing a subsequent AGC digitally is that you can measure the level (which sets the level of the overall waveform within the bandwidth at the detector) after more thorough filtering of the channel signal of interest and rejecting out of band interference (which could include other channels in a multiband receiver). This interference would reduce the signal of interest in band when AGC'd (power measured) earlier in the chain. Consider the above case earlier for the analog AGC where we have a very strong interference signal near the channel bandwidth of interest, so close that filtering in the analog would be costly. Assume too that we have a high dynamic range ADC with excess dynamic range allowing us to shift some of this filtering challenge to the digital. The analog AGC's purpose is to avoid excessive clipping at the ADC input (or other similar level sensitive components in the receiver chain), so with the strong interference our signal of interest will effectively be significantly attenuated, give the interference will set the full scale level at the ADC input.
So the Digital AGC's purpose is to further bring this signal back to the proper full scale level (and minimize subsequent receiver processing) by selecting the bandwidth of interest and removing the remaining interference digitally, and then once removed, bring the channel signal back to the proper full scale level (some determined back-off to minimize clipping in the waveform but also reduce the subsequent operations required assuming a fixed point design, meaning minimize datapath requirements).
Also we must be careful to consider all components between the adjusted signal level and the detector (assuming filtering prior to the detector means that the actual power in the signal affecting the components in line can be much higher leading to saturation, clipping and gain compression.

So analog AGC output is range of signal levels but digital AGC always
gives out a fixed level?

Not quite. Both analog and digital AGC provide a fixed level of channel power, as measured by the relevant power/signal detector used. What may more likely occur in the analog is strong interference controlling the AGC, so the AGC sets its level based on the interference and not our signal of interest. The AGC only measures power and does not recognize signal from interference: we help that with filter design.
There are also systems where the level is measured digitally and the front-end gain is controlled but it is important to monitor the total signal level at the input the ADC to optimize the trade between quantization noise and clipping, maximizing the SNR of the signal with the proper AGC set-point at the ADC input. (Some clipping is often actually better than backing off too far and introducing more quantization noise in the process).
The AGC loop consists of a detector and a way to adjust the power along with an appropriate loop filter as a control loop. The power adjustment can be done as gain or attenuation or both, depending on what is available in the design. The detection and gain adjust and loop filter can all be done digitally or analog or both and as suggested above it can be a multi-loop design, and the gain can also include stepped gain control.
For more info see these posts: How to adjust receiver gains to avoid saturation and quantization noise to optimise post digital processing?
Role of Power Metric and LPF in AGC(Automatic Gain Controller)
AGC noise level


Picking the correct filter for accelerometer data

As pointed out by @JohnRobertson in Bag of Tricks for Denoising Signals While Maintaining Sharp Transitions, Total Variaton (TV) denoising is another good alternative if your signal is piece-wise constant. This may be the case for the accelerometer data, if your signal keeps varying between different plateaux.

Below is a Matlab code that performs TV denoising in such a signal. The code is based on the paper An Augmented Lagrangian Method for Total Variation Video Restoration. The parameters $\mu$ and $\rho$ have to be adjusted according to the noise level and signal characteristics.

If $y$ is the noisy signal and $x$ is the signal to be estimated, the function to be minimized is $\mu\|{x-y}\|^2+\|{Dx}\|_1$, where $D$ is the finite differences operator.

function denoise()

f = [-1*ones(1000,1);3*ones(100,1);1*ones(500,1);-2*ones(800,1);0*ones(900,1)];
plot(f);
axis([1 length(f) -4 4]);
title('Original');
g = f + .25*randn(length(f),1);
figure;
plot(g,'r');
title('Noisy');
axis([1 length(f) -4 4]);
fc = denoisetv(g,.5);
figure;
plot(fc,'g');
title('De-noised');
axis([1 length(f) -4 4]);

function f = denoisetv(g,mu)
I = length(g);
u = zeros(I,1);
y = zeros(I,1);
rho = 10;

eigD = abs(fftn([-1;1],[I 1])).^2;
for k=1:100
    f = real(ifft(fft(mu*g+rho*Dt(u)-Dt(y))./(mu+rho*eigD)));
    v = D(f)+(1/rho)*y;
    u = max(abs(v)-1/rho,0).*sign(v);
    y = y - rho*(u-D(f));
end

function y = D(x)
y = [diff(x);x(1)-x(end)];

function y = Dt(x)
y = [x(end)-x(1);-diff(x)];

Results:


Designing a half band FIR filter with Scilab

In commerical filter design software, when designing an odd-tap half-band filter your center coefficient (the largest-valued coefficient) will have a value of 0.5. So if one of your filter's coefficient has a value of 0.000002, then that coefficient's value is one part in 250,000 compared to the largest coefficient. That super-small valued coefficient will have a negligible effect on the filter's freq response, so its value can be set to zero.

In MATLAB, try one of the following:

hn = remez(26, 2*[0 .23 .27 .5], [1 1 0 0], [1 1])
hn = firpm(26, 2*[0 .23 .27 .5], [1 1 0 0], [1 1])

(Notice the '2*' multiplier needed to comply with MATLAB's command syntax!) Examine the coefficients' values and see that none are exactly zero-valued. Now plot your filter's freq response (vertical axis in dB). Next, zero-out all coefficients whose absolute values are less than 0.001 using:

hn = hn.*(abs(hn)>0.001)

and plot the new freq response (in dB). Notice how the two freq responses are essentially identical.


Transform minimum phase FIR into linear phase FIR

It is generally impossible to transform a given minimum-phase FIR system into a linear phase FIR system with the same magnitude response. There is one special case for which this is possible, and that is if the zeros of the minimum phase system inside the unit circle haven even multiplicity. Because in that case you can, for each zero location, mirror half of the multiple zeros across the unit circle, which gives you a linear phase filter without changing the magnitude. But that's of course a quite artificial situation, so in general there's no way to obtain a linear phase filter from a minimum phase filter without changing the magnitude.

A way to transform a minimum-phase FIR system into a linear-phase FIR system with a squared magnitude is by simply doubling all zeros on the unit circle, and for each zero inside the unit circle adding a mirrored zero outside the unit circle. If $H(z)$ is a given $N^{th}$ order causal minimum-phase transfer function, the corresponding linear phase transfer function is given by

$$G(z)=H(z)H(1/z^*)z^{-N}\tag{1}$$

Note that the number of coefficients is (approximately) doubled and that

$$|G(e^{j\omega})|=|H(e^{j\omega})|^2\tag{2}$$


Relation between sawtooth Fourier coefficients and its DFT

Let me use the definitions of the Fourier series and the DFT as they are commonly used in signal processing. The Fourier series of a $T$-periodic function $f(t)$ is given by

$$f(t)=\sum_{k=-\infty}^{\infty}c_ke^{j2\pi kt/T}\tag{1}$$

where $c_k$ are the complex Fourier coefficients:

$$c_k=\frac{1}{T}\int_0^Tf(t)e^{-j2\pi kt/T}dt\tag{2}$$

The length-$N$ DFT of a sequence $x[n]$ is defined by

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{3}$$

and the inverse DFT is

$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{4}$$

If you sample $f(t)$ at $t=n\tau$ with $\tau=T/N$ you can approximate the integral $(2)$ by the following sum:

$$c_k\approx\frac{1}{T}\sum_{n=0}^{N-1}f(n\tau)e^{-j2\pi kn\tau /T}\tau\tag{5}$$

With $x[n]=f(n\tau)$ and $\tau=T/N$, Eq. $(5)$ can be rewritten as

$$c_k\approx\frac{1}{N}\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}=\frac{1}{N}X[k],\qquad |k|\le\left\lfloor \frac{N}{2}\right\rfloor\tag{6}$$

Equation $(6)$ says that the Fourier series coefficients $c_k$ of the $T$-periodic function $f(t)$ are approximated by the (scaled) DFT coefficients of the sequence $x[n]=f(n\tau)$, with $\tau=T/N$.

Of course this is also the case for your function. With

$$f(t)=A\left(\frac{t}{T}-\frac12\right),\qquad 0\le t< T\tag{7} $$

and $f(t)=f(t+T)$, we get

$$x[n]=f(nT/N)=A\left(\frac{n}{N}-\frac12\right),\qquad 0\le n< N \tag{8}$$

The Fourier coefficients of $f(t)$ are

$$c_k=\frac{jA}{2\pi k}\tag{9}$$

The DFT of $x[n]$ is given by

$$X[k]=\begin{cases}-\frac{A}{2},&\quad k=0\\
\frac{A}{e^{-j2\pi k/N}-1} = e^{j\pi k/N} \frac{jA/2}{\sin(\pi k/N)}=\frac{A}{2}(-1+j\cot(\pi k/N)),&\quad 0<k<N
\end{cases}\tag{10}$$

The first thing to notice is that $f(t)$ is an odd function (with a discontinuity at $t=0$), and, consequently, its Fourier coefficients are purely imaginary. On the other other hand, the sampled version of $f(t)$ is not odd because $x[0]\neq 0$; that's why we get a (relatively small) non-zero real part in the DFT coefficients. However, as predicted by Eq. $(6)$, the scaled DFT coefficients $X[k]/N$ still approximate the Fourier coefficients quite well. The zero real part of $c_k$ is "approximated" by the value $-A/2N$, which is small for large $N$. The figure below shows the imaginary part of $X[k]/N$ (red) and the (imaginary part of the) Fourier coefficients $c_k$ (blue):

The increasing difference between $c_k$ and $X[k]/N$ with increasing $k$ is due to aliasing caused by sampling $f(t)$. Note that $f(t)$ is not band-limited, so you will always have some amount of aliasing, no matter how fast you sample. The difference in the real part ($-A/2N$ vs $0$) is caused by the fact that $f(t)$ is odd, but $x[n]$ isn't because $x[0]=f(0^+)=A/2\neq 0$.

EDIT: Here comes the derivation of the DFT of $x[n]$ defined by $(8)$. From $(3)$ we have

$$X[k]=\frac{A}{N}\sum_{n=0}^{N-1}ne^{-j2\pi nk/N}-\frac{A}{2}\sum_{n=0}^{N-1}e^{-j2\pi nk/N}\tag{a}$$

Let's first consider the case $k=0$:

$$\begin{align}X[0]&=\frac{A}{N}\sum_{n=0}^{N-1}n-\frac{A}{2}\sum_{n=0}^{N-1}1\\&=
\frac{A}{N}\frac{(N-1)N}{2}-\frac{AN}{2}\\&=-\frac{A}{2}\end{align}\tag{b}$$

For $k\neq 0$, the second term in Eq. (a) equals zero. From the formula

$$\sum_{n=0}^{N-1}nq^n=\frac{(N-1)q^{N+1}-Nq^N+q}{(1-q)^2},\qquad q\neq 1\tag{c}$$

(which can be obtained by taking the derivative of the standard formula for a geometric series), and by noting that $e^{-j2\pi k}=1$, we get after some manipulations

$$\begin{align}X[k]&=\frac{A}{N}\sum_{n=0}^{N-1}ne^{-j2\pi nk/N}\\&=\ldots\\&=\frac{A}{e^{-j2\pi k/N}-1},\qquad k\neq 0\end{align}\tag{d}$$

from which all other representations shown in Eq. $(10)$ can be derived quite easily.


How to find system function, H(z) in the z-domain, given zero-pole plot of the system?

The $z^{-1}$ term in the $H(z)$ expression you provided is the expression of the zero at infinity:

$$\lim_{z\to\infty} \frac{1}{z} = 0$$

And in your specific case:

$$\lim_{z\to\infty} H(z) = G_{0}\frac{1}{z}\frac{(z-\frac{10}{9}j)(z+\frac{10}{9}j)}{(z-\frac{9}{10})(z+\frac{9}{10})} = 0$$

Generally, if the order of the denominator exceeds the order of the numerator, then you have 1 or more zeros at $\infty$


The textbook is right. A sine wave in the time domain has infinite energy since it continues over an infinite amount of time. When you transform into the Frequency domain all this energy is concentrated on a single (or two) frequency. Hence the value there is indeed infinite.

Of course these are all theoretical considerations. In the real world, ideal sine waves do not exist since they all have a beginning and an end.


Change of bin magnitude when phase of signal changes when doing FFT?

The FFT magnitude, for non-zero sinusoidal-like signals, will only be the same, given a change in phase relative to the window, for signals that are exactly periodic in the window width (the non-zero padded portion of the FFT aperture).


Distortion of a signal consisting of two periodic components

The question is trying to make you understand Aliasing

The low-pass filter is introduced (as a hint) to imply that the output is band-limited to half the Sampling frequency (Why ? See below for aliases).
This would allow you to understand how out of band signals can sneak into your sampled signals if left unchecked. Look at p5 here

You should read in detail about sampling.
Skipping the real math, as rule of thumb you can do as below to answer your questions quick

  a) for $f_0=1\,\text{kHz}$ the signal cannot be reconstructed

Once you sample your signal with a rate you create aliases at frequencies
$\pm f_{signal}\pm n*f_{sampling}$
where n can take any positive value

If you calculate all combinations for all your frequencies you end up with
[100 400 600 900 1100 1600 1900 2100 2900 ...]

Since you'd sampled with 1000Hz, and are using an (ideal) LPF at 500Hz, you'd take all components from the set above between [0-500Hz] in which case you end up with [100 400].

Do you see the 400Hz odd man out here?This is termed as the alias of the 600Hz signal.

So, for an input of [100 600] and you end up with [100 400].

  b) for $f_0=1.3\,\text{kHz}$ reconstruction works

Repeat steps above.

You'll find that the only components between [0-650Hz] are [100 600].

There you go !


Practical Gardner timing error implementation

If you have access to Mengali's book 'Synchronization Techniques for Digital Receivers', you will find the details how to choose it for a tracking loop, be it for frequency, phase or timing. Otherwise, just read about a standard discrete-time phase locked loop and see how the choice of paramters affect the loop bandwidth, which is what you are after.


WiFi Channelization for 40, 80 and 160 MHz

Your receiver does not know what BW you are using. Hence it will always be tuned to receive 20MHz BW (ch 36 or 40). So for the preamble part you need to transmit in channel bonding (2 x 20MHz) mode. Once the receiver has detected preamble it decodes the header to find out that data part is 40MHz. You need to turn on a bit in header for this (forgot the exact bit number) you can find it in standard. Same case for 80MHz BW. Preamble part is 4 x 20MHz.
Yes you need to continuously tune to 20MHz channels in the receiver. Otherwise no way of knowing what BW transmitter is sending. For fast synchronisation this is the feasible method. Once you detect preamble you need check for that header bit in both primary and secondary channels to check for 40MHz BW.


Concrete use of convolution algebraic properties

For part 3, there is an old trick dating from the time where computing was expensive. If you have a smoothed version $I_\sigma$ of a picture $I$, you can get a crude yet simple approximation of a  Laplacian from $I -I_\sigma$, illustrated here with the Gaussian:

It is indeed generalized in many unsharp filters, with other weights, or in multirate/multiscale (pyramid) operations.
For part 2, you can convolve twice by Gaussians with a single convolution of a single Gaussian. Yet the most classical (yet rarely understood) is that smoothing first then differentiating is the same as differentiating then smoothing, see for instance Laplacian of Gaussian operator. Many people believe that smoothing first reduces noise. However, one can combine derivatives and smoothing, and save a couple of operations, and even design innovative gradients:

For part 1, this is used a lot in proofs. This is something one uses implicitly very often, for instance in cases when you have an image that you need to convolve a lot of times (scale-space filtering for instance) vs a lot of images you convolve when a filter known a priori.
Related question:

How to choose the size of a Laplacian of Gaussian kernel for filtering images fast?


What is the rms value of a signal and how to use it

RMS is the root-mean-square value of a signal.  For a digitised signal, you can calculate it by squaring each value, finding the arithmetic mean of those squared values, and taking the square root of the result.

Loosely speaking, it represents the average "power" of a signal.  In electrical terms, all signals with the same RMS voltage, when applied across a resistor, will produce the same amount of heat.

Note that Matlab has an RMS function Y = rms(X) in the signal processing toolbox.  if you don't have that toolbox, it shouldn't be hard to write your own function.

Since the RMS value is only related to the signal amplitude, knowing it doesn't help you to create the waveform.  However, if you have a function that will generate the desired waveform, you could measure the RMS value of the result, and multiply by an appropriate scaling factor to give you the amplitude you want.

If "the noise contributes to 5% of the signal RMS value", then presumably the signal would only have an RMS value of 95% of the total before you added the noise.  But beware that RMS values don't simply add.  Adding noise of RMS amplitude 1 to a signal of amplitude 10 will not necessarily produce a new signal of RMS amplitude 11.


As this signal (function) exists inside an integration operator, and also assuming that there are no impulses residing at any discontinuities introduced by $u(t)$ step functions, then we can simply disregard the isolated discontinuity points and rely on the following:

$$[u(\tau)-u(\tau-4)][u(t-\tau)-u(t-4-\tau)] =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{for t < 0}}\\
&u(\tau)-u(\tau-t) &\scriptstyle{\text{for 0 < t < 4}}\\
&u(\tau-t+4)-u(\tau-4) &\scriptstyle{\text{for 4 <t < 8}}\\
&0 &\scriptstyle{\text{for t > 8}}\\
\end{align}
\end{cases}$$

Where the intervals are interpreted with the $t$ viewed as a parameter for an integrand argument function of a dummy variable $\tau$.

This convolution is reminiscent of convolving a rectangular pulse of duration 4 units with itself, whose output result should be a triangular pulse of duration 8 units, extending from $t=0$ to $t=8$ with an apex at $t=4$. Assuming that this reminiscence is correct, I would define the convolution result with this more apparent form as (defining the pulse as $p(t) = u(t) - u(t-4)$) :
$$y(t) = p(t)*p(t) =  \int_{-\infty}^{\infty} {p(\tau)p(t-\tau)d\tau} =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{for t < 0}}\\
&\int_{\tau=0}^{\tau=t}{1d\tau} = t &\scriptstyle{\text{for 0 < t < 4}}\\
&\int_{\tau=t-4}^{\tau=4}{1d\tau} = 8-t &\scriptstyle{\text{for 4 <t < 8}}\\
&0 &\scriptstyle{\text{for t > 8}}\\
\end{align}
\end{cases}$$

Where the limits of the integrals are deduced from the integrand arguments given at the first set of definition above: For example, for the interval $0<t<4$ the limits of integral are choosen based on the integrand $u(\tau) - u(\tau -t)$. Looking at this integrand one can see that for this to be nonzero (for a given fixed t) we need both $\tau > 0$, (from the first $u(\tau)$) and also $\tau < t$ from the second $u(\tau-t)$. So combining these two conditions requires that $0<\tau<t$ must be the limits of integration. The other cases are treated similarly.

On the other hand the operation is most easily understood and visualied by the help of an application of "flip & drag" technique to evaluate the intervals and to deduce the resulting integrand functions, instead of the way I presented as a piecewise definition form.


Clarification regarding discrete wavelet transform

Assuming that you have a sufficient number of samples, throwing away odd / even samples does not matter.

The DWT can be thought of as measuring time/frequency content with varying levels of time/frequency resolution.  For non frequency varying signals (like chirps), the odd and even samples both contain the same frequency and time content. (Unless we are talking about some explosive transient that occurs in one sample)  then you would need some a priori knowledge about that.

If the signal is of odd length, taking the odd samples would give you one more sample than the even samples but wouldn't make a difference in terms of information gained.


Converting a simple analog Butterworth filter into a digital filter

Converting the analog filter $H_a(s)$ into a digital filter $H_d(z)$ using the bilinear transform where T is the sampling period:

$\Large H_d(z) = H_a(s)\bigg|_{s=\frac{2}{T}\frac{z-1}{z+1}}$

Example:

Given a first order Butterworth filter

$H_a(s) = \frac{1}{1+RCs}$

$H_d(z) = H_a\bigg( \frac{2}{T} \frac{z-1}{z+1} \bigg)$

$H_d(z) = \frac{1}{1+RC\Big(\frac{2}{T}\frac{z-1}{z+1} \Big)}$

$H_d(z) = \frac{1}{1+RC(\frac{2}{T}\frac{z-1}{z+1})}$

$H_d(z) = \frac{1+z}{(1-2RC/T)+(1+2RC/T)z}$

$H_d(z) = \frac{1+z^{-1}}{(1+2RC/T)+(1-RC/T)z^{-1}}$

The coefficients of the denominator are the 'feed-backward' coefficients and the coefficients of the numerator are the 'feed-forward'  coefficients used to implement a real-time digital filter.


From your confusion of $x(t+a) \star h(t+a)$ vs $h(t) \star x(t+a)$
I guess that a little help on the argument manipulations on functions and convolutions could be appropriate here, moving in simple examples:

First let us express the usual simplistic case.  Consider the relation:
$$ y(t) = h(t) \cdot x(t) + g(t) \tag{1}$$

then manipulations on the argument $t$ is applied to all functions on both sides

$$ y(t+a) = h(t+a) \cdot x(t+a) + g(t+a)$$

or an arbitrary transform on $t$ would similary be:
$$ y(\phi(t)) = h(\phi(t)) \cdot x(\phi(t)) + g(\phi(t))$$

Now consider the case where two functions convolved to produce the third:
$$y(t) = \int_{-\infty}^{\infty} h(\tau) x(t-\tau) d\tau $$

which is abbreviated as
$$ y(t) = h(t) \star x(t) \tag{2} $$

Now be careful to interpret the case-2. The variable $t$ shows in all functions as an argument but you may not apply the transform on $t$ as you did with the case-1, so assume you have an arbitrary transform on $t$ as $\phi(t)$ then

$$ y(\phi(t)) \neq h(\phi(t)) \star x(\phi(t)) \tag{3} $$

For example, as in your case, if $\phi(t) = t+a$ then you get

$$ y(t+a) \neq h(t+a) \star x(t+a)$$
but
$$ y(t+a) = h(t) \star x(t+a) = h(t+a) \star x(t) $$

The justification of this can (only) be seen when you consider the integral definition of the convolution operator:

$$
\begin{align}
y(t) &= h(t) \star x(t) \\
& = \int_{-\infty}^{\infty} h(\tau) x(t-\tau)d\tau \\
y(t+a) & = \int_{-\infty}^{\infty} h(\tau) x((t+a)-\tau)d\tau \\
& = h(t) \star x(t+a) \\
\end{align}
$$

Note that since the live variable inside the integral only happens in one function ($x(t-\tau)$ in this case) then a change in $t$ will only affect one of them and you get:

$$ y(t+a) = h(t) \star x(t+a) $$
or from commutativity of convolution you get
$$ y(t+a) = h(t+a) \star x(t) $$

So this provides the answer you were looking for. However, it's not over. Because the following case represents an exception:

$$ y(-t) \neq h(t) \star x(-t) $$
but
$$ y(-t) = h(-t) \star x(-t) \tag{4} $$

So how to see this case-4. Again, using the integral definition :

Assuming that $y(t) = h(t) \star x(t)$, then compute the convolution between two new signals $g(t)=h(-t)$ and $z(t)=x(-t)$ as:
$$
\begin{align}
w(t) &= g(t) \star z(t) \\
& = \int_{-\infty}^{\infty} g(\tau) z(t-\tau)d\tau &g(\tau)=h(-\tau),z(t-\tau)=x(-(t-\tau)) \\
& = \int_{-\infty}^{\infty} h(-\tau) x(-(t-\tau))d\tau  &\text{ let } \tau'=-\tau \\
& = -\int_{\infty}^{-\infty} h(\tau') x(-(t+\tau'))d\tau'  &\text{ replace } \tau'  \text{ with } \tau \\
& = \int_{-\infty}^{\infty} h(\tau) x(-t-\tau) d\tau  \\
& = y(-t) \\
\end{align}
$$

hence we conclude that $h(-t) \star x(-t) = w(t) = y(-t) $. As stated before, you must always consult to the (explicit) integral definition to decide on the correct functions used in the convolution operator.


How do I recover the signal from an ECG image

Sorry, I use Mathematica, but it should be really easy to implement the idea into Matlab. I give the code anyway, so when my description is not detailed enough, you can get the rest from the code.

Basic idea is: you look at your image column-wise. Process every column of pixels separately. Note, that in the plot, I inverted the gray-values. So black is 1 and white is 0.
If you plot the (inverted) brightness-pixel-values you have basically only two situations. The first one is, when your column is not on a vertical grid line. The the plot looks like

The second situation is, where you are directly on a vertical grid line. Then the grid line influences the brightness of the whole column

But what you see is, that your dark EEG seems to be always the maximum. Therefore, the very complex algorithm is: Go through every column and take the position of the most black pixel.

img = ColorConvert[
   ImagePad[
    Import["http://i.stack.imgur.com/500Kg.jpg"], {{0, -20}, {0, 0}}],
    "Grayscale"];

Image[
 Transpose[Function[With[{m = Min[#]},
     Map[Function[{v}, If[v == m, 1, 0]], #]]] /@
   Transpose[ImageData[img, "Real"]]
  ]
 ]

Note, that I croped a bit of the right side of you image, since it was completely white there. The result is

Now you can join the points or interpolate them in every way you like and you get your EEG


Understanding the attenuation of a signal

Attenuation in this context refers to a reduction in amplitude of some signal component as it passes through a system (in your case, the system being the highpass filter. Linear filters are designed to pass certain frequency components while stopping (or attenuating) others. As the name implies, a highpass filter is designed to pass high-frequency content while attenuating low-frequency content.

As you noted, the tremor features that you're interested in are high in frequency and small in amplitude. Thus, they may be masked by low-frequency components in the original seismogram. A highpass filter can be used to accentuate the high-frequency tremor waves by attenuating any low-frequency waves that might also be present.


python scipy fft on numpy hanning window smears peaks

Is this meant to happen?

Yes, that is absolutely a well known effect of using any window function.
Taking a look through the Wikipedia article on window functions, we find that the rectangular window has the sharpest peak.
This is quite simply due to the fact that the rectangular window acquires the most data, meaning that it can distinguish between slightly different frequencies the best.

The problem with the rectangular window is that it has large side lobes.
In other words, if you measure a signal which doesn't exactly match one of the discrete Fourier transform (DFT) frequencies (a.k.a. bins), then with the rectangular window you get a lot of leakage into neighboring bins.
The formula for this is actually rather simple.
Suppose your sampling interval is $\delta t$ and you measure $N$ points.
Then the DFT frequencies are
$$f_n = n / N \delta t, \quad n\in[0,1,\ldots,N-1] \, .$$

Now suppose you measure a sinusoid with frequency $\xi$ so that the sampled values are $x_n = \cos(2 \pi \xi n \delta t)$.
The mod squared of the DFT values are then
$$ |X_k|^2 = \left( \frac{\sin(\pi(\xi - k))}{\sin(\pi(\xi-k)/N} \right)^2 \, . \tag{$*$}$$

We plot this function in Figure 1 for $\xi=24.4$.
Note that when $\xi$ is not an integer the DFT values are nonzero for all bins, and you can see from equation $(*)$ that the values of $|X_k|^2$ only fall off as one over the square of frequency.

Figure 1: DFT values for a cosine of frequency 24.4.
The blue line indicates the underlying function given in equation $(*)$ and the black dots are the values of that function evaluated at the integers.
In other words, the black dots are the DFT bin values.

The point of a window function is to make the leakage of power into bins near the real frequency (in this case 24.4) fall off faster, but this comes at the cost that the main peak is broadened.$^{[a]}$

  Does anyone else get this problem?

It's not a problem.
It's how windows are supposed to work.

  For my data, it seems a rectangular window is fine, but to me it seems that the Hanning window must be wrong here, as the integral of the entire spectrum should be the same (after I multiplied by 2 to correct amplitude), but in this case it seems to be different.

Applying a window changes the total power in your signal, as you have observed.
This is not a mistake at all.
It's pretty obvious that this has to happen because the integral (sum) of the power in the frequency domain has to equal the power you would have computed in the time doman (Parseval's theorem).
It's clear that multiplying e.g. a cosine by a window function changes the integral of the time domain signal.
This is called "loss" and must be taken into account when using windows!

$[a]$: I'm not actually showing this because I didn't bother to compute the points with the window, but the Wikipedia article shows many examples.


Do Spread Spectrum technologies combat time-selective fading of a channel?

I hadn't heard the term "time-selective fading" before, and after a bit of research it appears that it refers to the effects of doppler spreading.  Given that, the answer is "not really".

DSSS is typically extremely sensitive to carrier offset since the pseudo-random chipping sequence is usually quite long.  The long time means that a little carrier offset can lead to destructive interference.  This is dealt with by transmitting the data with a pilot.  It is easy to search for the pilot across multiple carrier offsets by correlating using FFT techniques.

FHSS does not really have any particular advantages or disadvantages that I am aware of when it comes to handling doppler.


Contours matching - finding contours displacement

Do you have to worry about a difference in scale between contours?  If not than you can simply find the centroid of each contour, and compute displacement by subtracting one from the other.  Then you can calculate the principal axes of the contours, and find the rotation angle between them.

If scaling is involved, then you can calculate the scale factor by taking the ratio of the corresponding principal axes.


Maybe you'd be better off using Bessel-Thompson filters. Maybe one could find an all-pass that approximately corrects a Butterworth filter, but it looks rather hard.


Designing filters vs. editing fft components

I'll here give another reason (perhaps the most important one) why implementing filters by manually changing frequency bins is usually a bad idea.
Note first that time domain convolution (simple FIR filter design) can also be implemented in the frequency domain and this procedure is known as FFT convolution ( despite its name the math operation here is multiplication and not convolution). The two techniques are equivalent but the latter is much trickier to implement correctly (although it's much faster than its time domain counterpart)
On the other hand, manually changing frequency bins (either by removing them entirely or just by reducing their amplitude),  means that you're discarding the original phase information, which in turn, can lead (and usually does) to inaccurate reconstruction of the modified time domain signal. Multiplying in the frequency domain does correspond to time domain convolution but, unlike true FFT convolution, we don't have a real filter kernel here so the results are less predictible.
Mathematically, perfect reconstruction is possible only when both original magnitude and phase information is retained. But since you're modifying both magnitude and phase, it's quite possible that there is no real signal (in time domain) that corresponds to your modified FFT data (at least not without significant digital artifacts being added to the reconstructed signal). That said, there are ways (e.g iterative procedures such as the Griffin-Lim algorithm ) that can be used to remedy this situation.


Short Time Fourier Transform analysis of newborn baby cry signal

Have a look at SAaC from LabROSA. You also might be interested in delta pitch without explicit pitch tracking.


The Fourier series of the cycloid can be expressed in terms of the Bessel functions of the first kind:
$$J_n(x)=\frac{1}{\pi}\int_0^{\pi}\cos(nt-x\sin t)dt,\qquad n\in\mathbb{Z}\tag{1}$$
Using the cycloid parameterization
$$y(t)=1-\cos t,\qquad x(t)=t-\sin t\tag{2}$$
which results in a period of $2\pi$ and a maximum value of $2$, the Fourier series of $y(t)$ as a function of $x$, referred to as $f(x)$, is given by
$$\bbox[#f8f1ea, 0.6em, border: 0.15em solid #fd8105]{
\begin{align}f(x)&=\frac32+\sum_{n=1}^{\infty}\frac{J_{n+1}(n)-J_{n-1}(n)}{n}\cos(nx)\\&=\frac32-2\sum_{n=1}^{\infty}\frac{J'_{n}(n)}{n}\cos(nx)\end{align}}\tag{3}$$
where $J'_n(x)$ is the derivative of $J_n(x)$ w.r.t. $x$.
The figure below shows the cycloid and its Fourier series approximation according to $(3)$ using the first $20$ coefficients in the sum:

Proof:
Let's consider the real-valued Fourier series. Since $f(x)$ is even, all sine coefficients vanish and we get
$$f(x)=a_0+\sum_{n=1}^{\infty}a_n\cos(nx)\tag{4}$$
with
$$a_0=\frac{1}{2\pi}\int_{-\pi}^{\pi}f(x)dx\tag{5}$$
and
$$a_n=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)\cos(nx)dx,\qquad n>0\tag{6}$$
With $dx=(1-\cos t)dt$ we obtain
$$a_0=\frac{1}{2\pi}\int_{-\pi}^{\pi}y(t)(1-\cos t)dt=\frac{1}{2\pi}\int_{-\pi}^{\pi}(1-\cos t)^2dt=\frac32\tag{7}$$
and
$$\begin{align}a_n&=\frac{1}{\pi}\int_{-\pi}^{\pi}(1-\cos t)^2\cos(nt-n\sin t) dt\\&=\frac{2}{\pi}\int_{0}^{\pi}(1-\cos t)^2\cos(nt-n\sin t) dt,\qquad n>0\tag{8}\end{align}$$
where the last equality follows from the fact that the integrand is even.
Expanding
$$(1-\cos t)^2=1-2\cos t+\cos^2 t=\frac32-2\cos t+\frac12\cos 2t\tag{9}$$
and using $\cos(\alpha)\cos(\beta)=\frac12[\cos(\alpha-\beta)+\cos(\alpha+\beta)]$, the integrand in $(8)$ can be rewritten as
$$\begin{align}(1-\cos t)^2\cos(nt-n\sin t)=\frac32\cos(nt-n\sin t) -&\\\big[\cos((n-1)t-n\sin t) + \cos((n+1)t-n\sin t)\big]+\\\frac14 \big[\cos((n-2)t-n\sin t) + \cos((n+2)t-n\sin t)\big]\tag{10}\end{align}$$
Plugging $(10)$ into $(8)$ and using the definition of the Bessel function $(1)$, we can write
$$a_n=3J_n(n)-2\big[J_{n-1}(n)+J_{n+1}(n)\big]+\frac12\big[J_{n-2}(n)+J_{n+2}(n)\big],\qquad n>0\tag{11}$$
This can be further simplified using the recurrence relation (10.6.1) for $J_n(x)$:
$$J_{n-1}(n)+J_{n+1}(n)=\frac{2n}{x}J_n(x)\tag{12}$$
which can be used to eliminate $J_{n-2}(n)$ and $J_{n+2}(n)$ in $(11)$, and which finally results in
$$a_n=\frac{J_{n+1}(n)-J_{n-1}(n)}{n},\qquad n>0\tag{13}$$
Since (10.6.1)
$$J_{n-1}(x)-J_{n+1}(x)=2J'_n(x)\tag{14}$$
where $J'_n(x)$ denotes the derivative of $J_n(x)$ with respect to $x$, the coefficients $a_n$ can also be expressed as
$$a_n=-\frac{2J'_n(n)}{n},\qquad n>0\tag{15}$$
Eqs $(7)$, $(13)$, and $(15)$ establish the result $(3)$.


What type of images are required for voxel classification?

In general it doesn't matter what the format is, because you will load it as array/matrix anyway.

That being said, I never saw preprocessing being done on DICOM images, usually they are transformed to nifty or something else before anything is being done to it. However, ignorance is not an argument.


Pretty much the same.
The intrinsic calibration tells you the relationship between stuff "out there" relative to the camera's frame of reference and the image the camera makes.  The extrinsic calibration tells you how to get from the camera's frame of reference to whatever your "global" frame of reference may be.


How does PSOLA pitch-shift work in real-time?

In PSOLA, pitch estimation is used to space out or overlap segments by integer multiples of the estimated pitch period in time.  Previous frames (or pitch period length sub-segments thereof) can be duplicated (as needed) in order to not leave any gaps between frames.  Segments can then be concatenated or cross-faded.  This can lengthen or shorted the composited result, and thus change the time duration.  Then this intermediate result can resampled to change the pitch (and possibly undo some or all of the time duration change), and to provide the proper number of samples needed for the output sample rate.

There will be some jitter, as frames are usually only shifted by a minimum of one whole pitch period, not by fractions thereof.  Buffering may be required to cover this time jitter, the amount related to the longest period (lowest pitch) expected to be handled.

The "clever math" might be in the choice and implementation of an accurate enough but low latency pitch estimation algorithm.  If the pitch cannot be estimated (due to consonants or polyphony, etc.), apparent quality may suffer.


It looks like the author implicitly makes the substitution $x(nT) \rightarrow x(n)$, but even if they were making this explicit, there are more clear ways to get to Eq. 4.3.

One way to do this is to denote the continuous (analog) and discrete time-domain signals as $x_a$ and $x$, respectively. Your text, using this convention, would then look like this (starting from Eq. 4.2):

$$
x_a(t) \xrightarrow{sampling} x_a(nT)
$$

and

$$
X(e^{j\omega})\ =\ \sum_{n = -\infty}^{\infty} x_a(nT) e^{-j\Omega nT}.
$$

(Note the right-hand-side exponent: it is $-j\Omega nT$, not $-j\Omega nt$ as it appears in your text; this looks like an error to me.) Then if you take into account the following equality (assuming uniform sampling):

$$
x(n)\ =\  x_a(nT), \quad -\infty < n < \infty,
$$

you end up with:

$$
X(e^{j\omega})\ =\ \sum_{n = -\infty}^{\infty} x(n) e^{-j(\Omega T)n},
$$

or (noting that $\omega = \Omega T$):

$$
X(e^{j\omega})\ =\ \sum_{n = -\infty}^{\infty} x(n) e^{-j\omega n}.
$$

By the way, this is a derivation from CFT to DTFT operation, not CFT to DFT one (as you say in the title of your question).


How do I calculate the step response of a discrete-time system?

Taking the $\mathcal Z$ transform of the impulse response,

$$H (z) = \dfrac{2}{1 - \frac 12 z^{-1}} - 1$$

Hence, the $\mathcal Z$ transform of the step response is

$$S (z) := \dfrac{H (z)}{1 - z^{-1}} = \cdots = \dfrac{3}{1 - z^{-1}} - \dfrac{2}{1 - \frac 12 z^{-1}}$$

Taking the inverse $\mathcal Z$ transform, we obtain the step response

$$s (n) = \begin{cases} 3 - 2 \left( \frac 12 \right)^n & \text{if } n \geq 0\\ 0 & \text{if } n < 0\end{cases}$$


Why histogram equalization fail in enhance images for bimodal histogram?

Yes it is true. A simple histogram equalization may fail for certain image intensity distributions. In your case there are some bins with excessive number of pixels. For other methods, try for example gamma and contrast corrections.


Looks like you need a general explanation of the discrete wavelet transform (DWT). DWT breaks a signal down into subbands distributed evenly in a logarithmic frequency scale, each subband sampled at a rate proportional to the frequencies in that band. The traditional Fourier transformation has no time domain resolution at all, or when done using many short windows on a longer data, equal resolution at all frequencies. The distribution of samples in the time and frequency domain by DWT is of form:

log f
 |XXXXXXXXXXXXXXXX  X = a sample
 |X X X X X X X X   f = frequency
 |X   X   X   X     t = time
 |X       X
 |X
 |X
  ----------------t

Single subband decomposition and reconstruction:

    -> high -> decimate -------------> dilute -> high
   |   pass    by 2      high subband  by 2      pass \
in |                                                   + out
   |                                                  /  =in
    -> low  -> decimate -------------> dilute -> low
       pass    by 2      low subband   by 2      pass

This creates two subbands from the input signal, both sampled at half
the original frequency. The filters approximate halfband finite impulse response (FIR) filters
and are determined by the choice of wavelet. Using Daubechies wavelets
(and most others), the data can be reconstructed to the exact original
even when the halfband filters are not perfect. Note that in the above scheme, the total amount of information (samples) stays the same throughout.

Decimation by 2: ABCDEFGHIJKLMNOPQR -> ACEGIKMOQ
  Dilution by 2: ACEGIKMOQ -> A0C0E0G0I0K0M0O0Q0

To get the logarithmic resolution in frequency, the low subband is
re-transformed, and again, the low subband from this transformation
gets the same treatment etc.

Decomposition:

    -> high -> decimate --------------------------------> subband0
   |   pass    by 2
in |                 -> high -> decimate ---------------> subband1
   |                |   pass    by 2
    -> low -> decim |                 -> high -> decim -> subband2
       pass   by 2  |                |   pass    by 2
                     -> low -> decim |
                        pass   by 2  |   .   down to what suffices
                                      -> .    or if periodic data,
                                         .     until short of data

Reconstruction:

subband0 -----------------------------------> dilute -> high
                                              by 2      pass \
subband1 ------------------> dilute -> high                   + out
                             by 2      pass \                /  =in
subband2 -> dilute -> high                   + dilute -> low
            by 2      pass \                /  by 2      pass
                            + dilute -> low
Start   .                  /  by 2      pass
here!   . -> dilute -> low
        .    by 2      pass

In a real-time application, the filters introduce delays, so you need
to compensate them by adding additional delays to less-delayed higher
bands, to get the summation work as intended.

For periodic signals or windowed operation, this problem doesn't exist -
a single subband transformation is a matrix multiplication, with wrapping
implemented in the matrix:

Decomposition:

|L0|   |C0  C1  C2  C3                | |I0|     L = lowpass output
|H0|   |C3 -C2  C1 -C0                | |I1|     H = highpass output
|L1|   |        C0  C1  C2  C3        | |I2|     I = input
|H1| = |        C3 -C2  C1 -C0        | |I3|     C = coefficients
|L2|   |                C0  C1  C2  C3| |I4|
|H2|   |                C3 -C2  C1 -C0| |I5|
|L3|   |C2  C3                  C0  C1| |I6|
|H3|   |C1 -C0                  C3 -C2| |I7|     Daubechies 4-coef:

     1+sqrt(3)        3+sqrt(3)        3-sqrt(3)        1-sqrt(3)
C0 = ---------   C1 = ---------   C2 = ---------   C3 = ---------
     4 sqrt(2)        4 sqrt(2)        4 sqrt(2)        4 sqrt(2)

Reconstruction:

|I0|   |C0  C3                  C2  C1| |L0|
|I1|   |C1 -C2                  C3 -C0| |H0|
|I2|   |C2  C1  C0  C3                | |L1|
|I3| = |C3 -C0  C1 -C2                | |H1|
|I4|   |        C2  C1  C0  C3        | |L2|
|I5|   |        C3 -C0  C1 -C2        | |H2|
|I6|   |                C2  C1  C0  C3| |L3|
|I7|   |                C3 -C0  C1 -C2| |H3|

C0, C1, C2, C3 are the "db2" lowpass FIR filter coefficients. Highpass
coefficients you get by reversing tap order and multiplying by
sequence 1,-1, 1,-1, ... Because these are orthogonal wavelets, the
analysis and reconstruction coefficients are the same.

A coefficient set convolved by its reverse is an ideal halfband lowpass
filter multiplied by a symmetric windowing function. This creates the
kind of symmetry in the frequency domain that enables aliasing-free
reconstruction. Daubechies wavelets are the minimum-phase, minimum
number of taps solutions for a number of vanishing moments (seven in
"db7" etc.), which determines their frequency selectivity.

I was asked to show the matrices for 6 coefficients, so here they are, made a bit larger for clarity but could be the same size as before too. Decomposition:

|L0|   |C0  C1  C2  C3  C4  C5                | |I0|
|H0|   |C5 -C4  C3 -C2  C1 -C0                | |I1|
|L1|   |        C0  C1  C2  C3  C4  C5        | |I2|
|H1|   |        C5 -C4  C3 -C2  C1 -C0        | |I3|
|L2| = |                C0  C1  C2  C3  C4  C5| |I4|
|H2|   |                C5 -C4  C3 -C2  C1 -C0| |I5|
|L3|   |C4  C5                  C0  C1  C2  C3| |I6|
|H3|   |C1 -C0                  C5 -C4  C3 -C2| |I7|
|L4|   |C2  C3  C4  C5                  C0  C1| |I8|
|H4|   |C3 -C2  C1 -C0                  C5 -C4| |I9|

Reconstruction:

|I0|   |C0  C5                  C4  C1  C2  C3| |L0|
|I1|   |C1 -C4                  C5 -C0  C3 -C2| |H0|
|I2|   |C2  C3  C0  C5                  C4  C1| |L1|
|I3|   |C3 -C2  C1 -C4                  C5 -C0| |H1|
|I4| = |C4  C1  C2  C3  C0  C5                | |L2|
|I5|   |C5 -C0  C3 -C2  C1 -C4                | |H2|
|I6|   |        C4  C1  C2  C3  C0  C5        | |L3|
|I7|   |        C5 -C0  C3 -C2  C1 -C4        | |H3|
|I8|   |                C4  C1  C2  C3  C0  C5| |L4|
|I9|   |                C5 -C0  C3 -C2  C1 -C4| |H4|

With:

C0 = 3.326705529500826159985115891390056300129233992450683597084705e-01
C1 = 8.068915093110925764944936040887134905192973949948236181650920e-01
C2 = 4.598775021184915700951519421476167208081101774314923066433867e-01
C3 = -1.350110200102545886963899066993744805622198452237811919756862e-01
C4 = -8.544127388202666169281916918177331153619763898808662976351748e-02
C5 = 3.522629188570953660274066471551002932775838791743161039893406e-02

More coefficient sets can be found here.


displaying speech data in log scale using Matlab

It doesn't really make sense to display an instantaneous AC waveform on a log scale - as demonstrated by your results. I googled for cooledit screenshots and found this.

To accomplish this you would want to keep your data linearly scaled and draw it normally. Then just label the axis at the desired dB points.


how to create a harmonic mask from fundamental frequency?

See this post which describes creating a notch filter using a delay and add:

Moving Average for Notch filtering

The filter used is generically given as

$$H(z) = 1+z^{-N}$$

This results in passing DC and placing zeros at the normalized radian frequency locations given by:

$$\omega_n = \frac{2\pi n}{N} + \frac{\pi}{N}$$

for all integers n and $\omega_n$ in the range of 0 to $2\pi$ to correspond to the primary frequency range of DC to the sampling rate.

This comes directly out of factoring the polynomial given by H(z), resulting in zeros on the complex unit circle at the angular frequencies given above. Because of the extra $\pi/N$ rotation on the unit circle, this filter works well for filtering odd harmonics based on prudent choice of N and therefore the sampling rate as an integer multiple of the notch location.

If instead a delay and subtraction is used:

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

The result would be a zero (notch) at DC and and placing zeros at the normalized radian frequency locations given by:

$$\omega_n = \frac{2\pi n}{N}$$

Which comes from factoring $H(z) = 1-z^{-N}$ which results in the Nth roots of unity (zeros evenly spaced around the unit circle). This approach conveniently allows for placing a notch at every harmonic of a selected frequency.  To introduce the more advanced approaches given below, it is useful at this point to visually plot the zeros in the z domain for this case, using the example $H(z) = 1- z^{-8}$, along with the associated frequency response:

Both of these approaches are referred to as "Comb Filters" and notice in this case how the notch occurs at a fundamental frequency and every integer multiple of that frequency. The requirement here is the sampling rate be an integer multiple of the notch frequency. The filter has linear phase which is important if phase distortion is of concern, but a potential drawback is the gradual loss in frequency components at other than the distinct harmonics.

An improved approach for the amplitude response would be to extend the concept of the second order notch filter which I have derived at this link given below.

Transfer function of second order notch filter

This will result in much flatter amplitude response over the pass frequencies, but the tight notches come at the cost of increased phase distortion in proximity of the notches which may be an issue depending on the application.

To extend this to multiple notches at integer multiples note that the background of the 2nd order notch filter is the simpler DC notch filter given by

$$ H(z) = \frac{1+a}{2}\frac{(z-1)}{(z-a)} $$

Which has a zero at $z=1$ and a pole placed inside the unit circle at $z = a$ in close proximity to that zero; the closer the placement the tighter the bandwidth but also the more precision (bit size) is required to achieve the desired performance and maintain stability. Notice below the comparison to a single zero as done in the Comb Filter implementations (with $N=1$) to an example case of the DC notch filter ($a=0.95$ for the example shown):

At this point you may see where I am going with this. We can combine the concept of the DC Nulling filter with the comb filter to achieve tight notches at integer harmonic locations by placing poles in close proximity to the zeros given in the Comb Filter example as depicted in the graphic below.

Following similar polynomial factoring relationships into the complex plane, this is realized with the following transfer function below:

$$H(z) = K\frac{1-z^{-N}}{1-a^Nz^{-N}}$$

The amplitude normalization factor K if desired is determined by solving for the value in the passpand $H(z=\omega_n)$ where $\omega_n = \frac{2\pi n}{N} + \frac{\pi}{N}$ and inverting this, which ends up being:

$$K = \frac{1+a^8}{2}$$

The result for N= 8 and a = 0.95 is shown in the plot below.

The possible challenges with this approach is the significant passband phase distortion as N increases, and the amount of precision required in implementation.

If the phase distortion was an issue and the notches need not be as tight as the previous IIR implementation can provide, then a third option is to combine the linear phase notch response given by $H(z) = 1-z^{-N}$ cascaded with a simple 3-tap linear phase compensator similar to what I detailed in this link (the compensator shown is actually a raised minus cosine response; a raised positive cosine response would be the approach for this case, following the same process outlined at that post.)

how to make CIC compensation filter

In particular see the plot in that post for the interpolated CIC compensator which is implmented by simple zero insert (meaning to interpolate by 8, simply insert 7 zeros between each filter coefficient in the 3 tap compensator). This won't achieve the flatness and sharp nulls as done with the IIR approach above but will have significantly improved flatness over a wider passband range than the simple notch filter approach AND with linear phase. The cascade of the simple delay and subtract comb filter with the interpolated 3 tap FIR compensator (given by the two responses shown below) results in a significantly flatter passband response with minimal added complexity.


fundamental question regarding bandwidth

A single frequency is a pure sinusoid, and thus unmodulated (no BPSK).  Once you add any modulation, you no longer have a pure sinusoid.  The FT of this distorted sinusoid will have a non-zero bandwidth from the sidebands required to change the sinusoid from just pure to somewhat distorted (same as modulated).


How to get RMS to Frequency chart for discrete samples of acceleration?

If your time domain data is in units of vertical acceleration, then you can simply multiply the normalized (divide by N) raw fft values times the complex conjugate of itself and then take the square root. Each result will be the rms acceleration at that frequency.

Note, if your data is sampled at 100 Hz, you are going to see the content from 0 to 50 Hz of all frequencies that are within the analog bandwidth of your sampler. Any higher analog frequencies will "alias" into the 0 to 50 Hz band. I am not sure what you have for an anti-alias filter, but it certainly can't be a "brick-wall" filter right at 50 Hz. A properly designed filter may have hope of sufficiently rejecting signals above 60 Hz, in which case you can trust your data from 0 to 40 Hz (for example; I am more illustrating the need to understand this before interpreting your data).

Regarding "Normalized FFT"; some FFT implementations will have this built in while others you will need to normalize. You can test this by taking an FFT of 512 ones [1 1 1 1 ...1] and see if the first sample of the FFT is 512 or 1. If 512 then you need to divide the FFT result by the number of samples (512) to normalize the result.  And by "raw FFT" data I am referring to the 512 complex FFT samples, not the FFT frequency spectrum you plotted above which I assume is the power spectrum which already performed the complex conjugate multiplication and is only showing the positive frequency axis, likely with a dB vertical scale.

Usually these limits versus frequency are also normalized to be per Hz bandwidth. The FFT noise bandwidth is 1 bin (without further windowing which is another topic of importance if you start seeing dominant tones). So if you sample at 100 Hz with 512 bins, then each FFT bin has a noise bandwidth of $100/512= 0.1953$ Hz.  A more accurate measurement would average the power (the complex conjugate multiplication of the normalized FFT) over 5.12 bins before taking the square-root to get the average rms value per Hz bandwidth. A simpler and more conservative approach (will bound the actual max limit) would be to multiply the power in each bin by 5.12 before taking the square-root.

You may see spikes above the threshold in the time domain, but the limit is rms and over frequency, so those results will expect to be less.


Phasor notation of a complex number with apparently non-constant radius

Hint: a phasor is defined by its magnitude and phase. The magnitude can be obtained by computing

$$
    \sqrt{\left(\mathcal{Re}\left\{x_1[n]\right\}\right)^2
      + \left(\mathcal{Im}\left\{x_1[n]\right\}\right)^2}
$$

And the phase can be obtained with

$$
    \tan^{-1}{\frac{\mathcal{Im}\left\{x_1[n]\right\}}{\mathcal{Re}\left\{x_1[n]\right\}}}
$$


Confusion about speech being piecewise stationary

In this context, stationary usually means that the spectral density is (theoretically) the same throughout the observation window.
So, piecewise stationary means that if you divide the signal into short enough time frames, the signal in each frame can be described, accurately enough, as sum of sinusoids in presence of white noise (which is stationary as well).


Where to get the Kodak PhotoCD raw image dataset?

The dataset appears to be in the public domain, and has been posted at this site here.

The pictures (24 in total) were originally released in 1991 on a CD.  The original media with the actual image data could be requested via interlibrary loan; here is a WorldCat link to the repositories holding this item. The CD was originally distributed by Apple Computer.  The raw version of the files do not appear to be available online.


ifft of signal without phase information

This is tricky. The phase information impacts the time domain form very much. For example white noise and a delta impulse have exactly the same amplitude spectrum, the only difference is in the phase.

Without any further assumptions or extra knowledge this can't be done. For example if you know that the signal is stationary, you could just try a random phase. If it's the impulse response of a causal system, you could try a minimum phase.


$$ \mathscr{F} \Big\{ x(t) \Big\} \triangleq X(f) \triangleq \int\limits_{-\infty}^{\infty} x(t) \, e^{-j2 \pi f t} \ \mathrm{d}t $$

$$ \mathscr{F}^{-1} \Big\{ X(f) \Big\} \triangleq x(t) = \int\limits_{-\infty}^{\infty} X(f) \, e^{j2 \pi f t} \ \mathrm{d}f $$

Assuming $x(t)$ is real (which means that $R_x(\tau)$ is also real).

$$\begin{align}
 S_x(f) &\triangleq \mathscr{F} \Big\{ R_x(\tau) \Big\} \\
\\
 &= \int\limits_{-\infty}^{\infty} R_x(\tau) \, e^{-j2 \pi f \tau} \ \mathrm{d}\tau \\
\\
 &= \int\limits_{-\infty}^{\infty} \left( \int_\limits{-\infty}^{\infty}x(t)x(t-\tau) \mathrm{d}t \right) \ e^{-j2 \pi f \tau} \ \mathrm{d}\tau \\
\\
 &= \int\limits_{-\infty}^{\infty} \int_\limits{-\infty}^{\infty}x(t)x(t-\tau) \ e^{-j2 \pi f \tau} \,\mathrm{d}t \,\mathrm{d}\tau \\
\\
 &= \int\limits_{-\infty}^{\infty} \int_\limits{-\infty}^{\infty}x(t)x(t-\tau) \ e^{-j2 \pi f \tau} \,\mathrm{d}\tau \,\mathrm{d}t \\
\\
 &= \int\limits_{-\infty}^{\infty} x(t) \int_\limits{-\infty}^{\infty}x(t-\tau) \ e^{-j2 \pi f \tau} \,\mathrm{d}\tau \,\mathrm{d}t \\
\\
 &= \int\limits_{-\infty}^{\infty} x(t) e^{-j2 \pi f t} \int_\limits{-\infty}^{\infty}x(t-\tau) \ e^{-j2 \pi f (\tau-t)} \,\mathrm{d}\tau \,\mathrm{d}t \qquad \text{let } u \triangleq \tau-t\\
\\
 &= \int\limits_{-\infty}^{\infty} x(t) e^{-j2 \pi f t} \int_\limits{-\infty}^{\infty}x(-u) \ e^{-j2 \pi f u} \,\mathrm{d}u \,\mathrm{d}t \\
\\
 &= \int\limits_{-\infty}^{\infty} x(t) e^{-j2 \pi f t}  \,\mathrm{d}t \int_\limits{-\infty}^{\infty}x(-u) \ e^{-j2 \pi f u} \,\mathrm{d}u \\
\\
 &= \mathscr{F} \Big\{ x(t) \Big\} \cdot \mathscr{F} \Big\{ x(-t) \Big\} \\
\\
 &= X(f) \cdot X(-f) \\
\\
 &= X(f) \cdot \big( X(f) \big)^*  \\
\\
 &= \Big|X(f)\Big|^2  \\
\end{align}$$


For me, only the case of an even number of members of the sum is clear. How to prove this in general terms?

ah, "sum orthogonality proof" is something you'll want to google. You can do it (under a relaxed meaning of "proof") very geometrically. Understand each part in your sum as vector in a plane:

take the first vector: $e^{0}\equiv 1$. It goes from point 0+0j to 1+0j.
attach the second vector to the tip of that vector. Note the inner angle between these two. Also note how the position of the tip is the partial sum of the first two elements from your sum.
attach the third to the tip of the second. Note that the angle between these two is the same.
Repeat until you're out of vectors.

You clearly followed a polygonal chain, where each of the $N$ edges had length 1. Now, what is the condition on the angles for such a polygonal chain to be a closed polygon?


How do you properly organize data to compute multiple (independent) recursive filters at the same time taking advantage of SIMD instructions?

usually, SIMD instructions work by loading the values you want to apply the same operation to into a wide register, and then applying the vectorized instruction to it.

Now, collating values from all over memory to end up in one say 256-bit wide register, just to do 8 multiplications at once, probably doesn't pay.

If you, however, are able to decide the in-memory layout of these parallel streams $x_i$ to be interleaved, i.e.

\begin{matrix}
x_0[0] & x_1[0] & x_2[0] &\cdots&x_7[0]&x_0[1] & x_1[1]&\cdots &x_7[1] & x_8[2]&\cdots
\end{matrix}

then you could load the consecutive values into one SIMD register without need for any memory or register shuffling; that would be desirable.

Again, in your specific case I wouldn't expect too much from optimization for SIMD: If a single-tap IIR is already the bottleneck, your application is not CPU, but memory bandwidth bound.

If you just think you can speed up your application significantly by speeding up this especially cheap operation: I'm afraid that's rather unlikely; if it's trivial to parallelize on an instruction level, your compiler stands a good chance of doing exactly that, and really, one multiply-accumulate per output sample is pretty idyllic, so I doubt there's much to be done, considering the cost of rearranging numbers to make that work (your API defines samples to not be in interleaved order, and your compiler can't do anything about that).

There's other aspects to be looked at; among these are questions like whether double is really necessary for what is presumably fixed-point numbers with less than 32 bit of dynamic range.


proof of Autocorrelation property of DFT

You did not consider the conjugate property of the DFT.

The complex conjugate property of DFT states that
$\rm{DFT}\{x^{*}[n]\} = X^{*}[-k]$.
This implies that exponent in the second line of your proof should be $-k$. Which then allows for the DFT of $x[n]$ as $X[k]$.
A number multiplied by its conjugate give the square of its absolute value.

And DFT is discrete in the transform domain, hence its parameter should be in brackets rather than parentheses.


Why is the Fourier transform of a Dirac comb a Dirac comb?

I believe that the fallacy is to believe that a Dirac comb is localized in time. It isn't because it is a periodic function and as such it can only have frequency components at multiples of its fundamental frequency, i.e. at discrete frequency points. It can't have a continuous spectrum, otherwise it wouldn't be periodic in time. Just like any other periodic function, a Dirac comb can be represented by a Fourier series, i.e. as an infinite sum of complex exponentials. Each complex exponential corresponds to a Dirac impulse in the frequency domain at a different frequency. Summing these Dirac impulses gives a Dirac comb in the frequency domain.


Detection of a Circle in Noisy Image Data

Here is my solution, it is close to @Yoda's idea, but I changed some steps.

Mark all pixels such that there are at least 6 pixels in their 7x7 neighborhood
Remove all blobs, but the largest
Fill holes
Apply edge detection
Find circle using Hough transform

Here is the relevant Matlab code. I am using Hough transform for circles .m file in my code.

function FindCircle()
    close all;
    im = imread('C:\circle.png');
    im = im(:,:,2);

    ims = conv2(double(im), ones(7,7),'same');
    imbw = ims>6;
    figure;imshow(imbw);title('All pixels that there are at least 6 white pixels in their hood');

    props = regionprops(imbw,'Area','PixelIdxList','MajorAxisLength','MinorAxisLength');
    [~,indexOfMax] = max([props.Area]);
    approximateRadius =  props(indexOfMax).MajorAxisLength/2;

    largestBlobIndexes  = props(indexOfMax).PixelIdxList;
    bw = false(size(im));
    bw(largestBlobIndexes) = 1;
    bw = imfill(bw,'holes');
    figure;imshow(bw);title('Leaving only largest blob and filling holes');
    figure;imshow(edge(bw));title('Edge detection');

    radiuses = round ( (approximateRadius-5):0.5:(approximateRadius+5) );
    h = circle_hough(edge(bw), radiuses,'same');
    [~,maxIndex] = max(h(:));
    [i,j,k] = ind2sub(size(h), maxIndex);
    radius = radiuses(k);
    center.x = j;
    center.y = i;

    figure;imshow(edge(bw));imellipse(gca,[center.x-radius  center.y-radius 2*radius 2*radius]);
    title('Final solution (Shown on edge image)');

    figure;imshow(im);imellipse(gca,[center.x-radius  center.y-radius 2*radius 2*radius]);
    title('Final solution (Shown on initial image)');

end


What you're probably looking for is the inverse (discrete Fourier) transform of $|X(e^{j\omega})|^2$, i.e., of $|X(z)|^2$ for $|z|=1$. Otherwise, as pointed out in Fat32's answer, $|X(z)|^2$ is generally no valid $\mathcal{Z}$-transform.

Note that a valid $\mathcal{Z}$-transform must have the form of a power series:

$$X(z)=\sum_{n=-\infty}^{\infty}x[n]z^{-n}\tag{1}$$

Take as a simple example the $\mathcal{Z}$-transform of $x[n]=\delta[n]+\delta[n-1]$:

$$X(z)=1+z^{-1}\tag{2}$$

Then $|X(z)|^2$ is given by

$$|X(z)|^2=\left|1+z^{-1}\right|^2=1+2\text{Re}\{z^{-1}\}+\left|z^{-1}\right|^2\tag{3}$$

which is not of the form $(1)$, hence it is no valid $\mathcal{Z}$-transform.

Since $|X(e^{j\omega})|^2=X(e^{j\omega})X^*(e^{j\omega})$, you need to find the IDFT of $X^*(e^{j\omega})$:

$$\frac{1}{2\pi}\int_{-\pi}^{\pi}X^*(e^{j\omega})e^{jn\omega}d\omega=\frac{1}{2\pi}\left[\int_{-\pi}^{\pi}X(e^{j\omega})e^{-jn\omega}d\omega\right]^*=x^*[-n]\tag{4}$$

Consequently, the IDFT of $|X(e^{j\omega})|^2$ is given by $x[n]\star x^*[-n]$, where $\star$ denotes convolution.

It's also good to know what this means in terms of the $\mathcal{Z}$-transform. The $\mathcal{Z}$-transform of $x^*[-n]$ can be derived as follows:

$$\sum_nx^*[-n]z^{-n}=\sum_nx^*[n]\left(\frac{1}{z}\right)^{-n}=\left[\sum_nx[n]\left(\frac{1}{z^*}\right)^{-n}\right]^*=X^*\left(\frac{1}{z^*}\right)\tag{5}$$

So the $\mathcal{Z}$-transform equivalent of $|X(e^{j\omega})|^2$ is $X(z)X^*(1/z^*)$, as they correspond to the same time domain sequence. As shown above, the function $X(z)X^*(1/z^*)$ is a valid $\mathcal{Z}$-transform, unlike $|X(z)|^2$.


What is the name of the following 1st order high-pass IIR filter?

This is a first-order discrete-time Butterworth filter, which can be obtained via the bilinear transform from the following analog prototype filter (with normalized cut-off frequency):

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

The bilinear transform replaces $s$ by $k\frac{z-1}{z+1}$, where $k$ is chosen such that the discrete-time filter has the specified cut-off frequency $f_c$, resulting in

$$k=\frac{1}{\tan(\pi f_c T)}\tag{2}$$

where $T=1/f_s$ is the sampling interval.

So if you replace $s$ in $(1)$ by

$$s=\frac{1}{\tan(\pi f_c T)}\frac{z-1}{z+1}\tag{3}$$

you obtain the corresponding discrete-time filter with the coefficients given in your question.

EDIT:

OK, here is how you obtain your final coefficients. Let $a=\tan(\pi f_c T)$. Plugging $(3)$ into $(1)$ then gives for the transfer function of the discrete-time filter

$$\begin{align}G(z)=\frac{\frac{1}{a}\frac{z-1}{z+1}}{1+\frac{1}{a}\frac{z-1}{z+1}}&=\frac{z-1}{a(z+1)+z-1}\\&=\frac{z-1}{(1+a)z-(1-a)}\\&=\frac{1}{1+a}\frac{z-1}{z-\frac{1-a}{1+a}}\end{align}\tag{4}$$

which is exactly in the given form with $c_1=1/(1+a)$ and $c_2=(1-a)/(1+a)$.


Which window function to use for QPSK signal?

For purpose of recovering the clock, no window function is needed. Also to clarify, your signal is downsampled to 10GS/s which is samples per second, not "Symbols per second" as you specified the symbol rate to be 2.4 GBuad which is 2.4 G Symbols/second.

For QPSK clock recovery, I recommend considering the Gardner Loop in which case you would filter and downsample to 2 samples per symbol (4.8 GS/s). You can see the link below for more information on the Gardner Loop implementation and considerations.  I have seen very efficient implementations for closing the timing loop using polyphase filters; let me know if you would like more information on that (and update your title to not mention windowing) and I can add that here.

You will likely also have a carrier offset that you will need to track and remove. The Gardner is fairly insensitive to carrier offset, so once timing is acquired you can downsample to one sample per symbol (at the best decision locations) and use those complex samples to determine the carrier offset and correct for that in a carrier tracking loop (such as a digital Costas Loop).

See: Gardner Timing Recovery for Repeated Sybmols


How does cascading biquad sections for higher order filters work?

It's the last thing you said ("Or does the output of the first filter feed as x_in in to the second filter and so on?"). The idea is simple: you treat the biquads as separate second-order filters that are in cascade. The output from the first filter is the input to the second, and so on, so the delay lines are spread out among the filters. If you need to optimize the structure in a memory-constrained environment, you can note that adjacent biquads have redundant delay memory (i.e. the last few output samples of stage 1 are the same as the last few input samples of stage 2, so you don't have to store them separately like you would if you just implement the filters in isolation).


What are some suitable methods for removing low-frequency line noise from measured data prior to system identification/deconvolution?

Taking a quick look at your acquisition: this doesn't look too terrible in terms of noise, however, the acquisition window looks way too short. With a 60ms window your frequency resolution is only 160 Hz. From the looks of it, you will also get a massive truncation error. As a rough rule of thumb, you want the acquisition window long enough so that the impulse response visually disappears in the noise.

Specific answers (although they may not help given the truncation error)

If it's really line (like in AC supply voltage hum) noise you can measure the exact location and notch it out for the fundamental and the harmonics.
Sure. Take the mean and subtract. However, again, the window is way too short to get decent low frequency resolution.

A bunch of tips

Make the acquisition window long enough to capture the complete impulse response
Take an acquisition with no excitation. That's your noise floor
Take an acquisition with excitation. That's your signal
Compare signal spectrum to noise floor spectrum. Calculate SNR as a function of frequency
You want at least 10 dB SNR for a decent measurement. Frequencies must either meet that or be so low in level that they don't matter (< at least 20 dB of "the bulk")

If that can't be done there are more esoteric methods like Wiener filtering and/or reconstructive modelling, but it's certainly preferable to clean up the measurement process.


How to detect the singular value with signal processing technology?


Multiple biquads vs higher order filtering

A cascade of 2nd order biquadratic (biquad) sections is a commonly used approach for  implementing higher order FIR/IIR filters in fixed point systems (Once you determine that you need an IIR filter, this would be a topology to consider). A big consideration for the biquad structure is stability: higher order IIR filters are prone to instability in fixed point systems due to coefficient rounding combined with long feedback delay paths (causing poles to go outside of the unit circle). The absence of any feedback between the individual cascaded sections (decoupling the poles and zeros between sections) results in a more stable structure that will be less sensitive to coefficient quantization.
I do not know of other advantages, but the numerical robustness is a very serious consideration especially  in higher order structures (and fixed point systems).


Fourier Transform Signals - Time Transformations

HINT:

Note that $x(t)+x(t-1)$ has the required shape. Now it's only a matter of shifting, time scaling, and amplitude scaling to arrive at the exact expression for $z(t)$.


Should you include division by coefficient of variation in WV reconstruction for its normed 1D signal?

There is no benefit of dividing the new improved signal by the covariance of the original ECG and its improved signal like the following code; by the two decimal accuracy in Corrgrams and other more precise tools (eigen cor)

% Covariance http://dsp.stackexchange.com/q/35291/5633
% https://en.wikipedia.org/wiki/Coefficient_of_variation
utp = utp / sum( cov(u,utp) )';

Fig. 1 Left-hand-side figure without division by covariance and Right-hand-side figure with division by covariance - about marginal dependence and independence of the cases

Other inspected methods where I see no difference

Eigen cor matirx of the two
-M^{-1} about conditional independence and dependence

Please, let me know how to get more decimal digits in R corrgram panel.cor; I am trying to replace the problem by Didzis' p-value corrgram discussed there but having problems with finding the upper bound for the amount of points to maintain significant p-value.


How is the complexity of applying Short term fourier transform over n samples with a moving window of size m amount to O(nmlogm)

evidently the author is performing an FFT $n$ times (or, if the window hop is not 1 sample, it would be proportional to $n$) and each FFT costs something proportional to $m \log(m)$).  all this is quite normal.


The sequence produced by the Voss algorithm is a sequence of random numbers whose power spectral density decays with frequency as $1/f$. In other words, if you find the magnitude spectrum of the sequence (in other words, the magnitude of the DFT of the sequence), you'll see that it decays roughtly as $1/f$.

As to the meaning of the sequence: it doesn't have any intrisic meaning. One way to think about it is as noise with a particular frequency spectrum. In other words, if you normalize the sequence and feed it to your sound card, you'll hear a certain kind of noise. This noise can be added to other signals, such as audio, to give it specific properties.

If all you need is $1/$ noise, there are other methods to generate it; see Pink ($1/f$) pseudo-random noise generation.


Based on your exploration of the topics I can make the following explanation, mostly following the approach in Statistical Digital Signal Processing and Modeling, HAYES.

Linear Time Invariant (LTI) filtering is associated with the following LCCDE (Linear Constant Coefficient Difference Equation) with zero initial conditions and assuming a causal solution; i.e., interested in the solution for $n \geq 0$ in general.

$$ \sum_{k=0}^{N} a_k y[n-k] = \sum_{k=0}^{M} b_k x[n-k]  \tag{1} $$

where $a_k$ and $b_k$ are the coefficients of the LCCDE (mostly $a_0 = 1$), and $x[n]$ is the input, while $y[n]$ is the output of the filter (aka solution of LCCDE).

The filter is associated with an impulse response $h[n]$ (coefficients of the LTI filter) which is the solution of Eq(1) when $x[n] = \delta[n]$; $x[0] = 1$ and $x[n] = 0$ for all $n> 0$.

$$ x[n] \longrightarrow \boxed{ h[n] } \longrightarrow y[n] $$

Impulse response can be used to compute the solution $y[n]$ of the LCCDE in Eq(1) for a given input $x[n]$, through what's called a convolution sum :

$$ y[n] = h[n] \star x[n] = \sum_{k=-\infty}^{\infty} h[k] x[n-k] \tag{2} $$

Z-transform of the impulse response is the following and is useful as an algebraic transform representation of the LTI filter system.

$$ H(z) = \frac{ \sum_{k=0}^{M} b_k z^{-k} }{\sum_{k=0}^{M} a_k z^{-k}} \tag{3}$$

From which the output $y[n]$ can also be expressed as
$$ Y(z) = H(z) X(z) \tag{4} $$

which is nothing but Z-transform applied to Eq(2).

ARMA-(p,q) modeling :

Consider an LTI filter with LCCDE coefficients $a_k$ and $b_k$, as given in Eq(1), then if the input of this filter is a white-noise $v[n]$ with unit variance, then the output of this filter is called (by definition) Auto-Regressive Moving Average process with orders (p,q), abb as ARMA-(p,q), the associated LCCDE is written like

$$ \sum_{k=0}^{p} a_k x[n-k] = \sum_{k=0}^{q} b_k v[n-k]  \tag{5} $$

Note that with the filtering convention in Eq(1), $x[n]$ was the input and $y[n]$ was the output, but withing the ARMA-(p,q) convention in Eq(5), $v[n]$ is the input and $x[n]$ (the ARMA process) is the output. The AR order is $p$ (associated with the lefthand side, recursive part, of the LCCDE) and the MA order is $q$ (associated with the right hand side of the LCCDE).

Now given some WSS RP $x[n]$, then you can model this process as an ARMA-(p,q) by finding the set of coefficents $a_k$ and $b_k$ according to Eq(5), which would minimize the mean square error between the modeling ARMA-(p,q) filter output $\hat{x}[n]$ and the given process $x[n]$ to be modeled. The solution for $a_k$ are generated by Modified Yule-Walker Equations and $b_k$ are found in a more complex manner following MA modeling techniques etc.

Linear Prediction :

Given a WSS random process $x[n]$, we want to estimate its current value $x[n]$ by  using its $p$ past values  $x[n-1],x[n-2],...,x[n-p]$ via a linear combination :

$$ \hat{x}[n] = \sum_{k=1}^{p} c_k x[n-k] $$

I have used these $c_k$ as the linear prediction coefficients; it can be shown that these coefficients are related to all-pole signal model Prony coefficients, or an ARMA-(p,0) model coefficients as $c_k = -a_k$.

Rewriting Eq(5) with $q=0$,and $a_0 = 1$, as

$$ x[n] = - \sum_{k=1}^p a_k x[n-k] + b_0 v[n] $$

with a suitable $b_0$, reveals that, if the predicted process $x[n]$ was generated via an ARMA-(p,0) model using Eq(5), then the linear prediction coefficients $c_k$ which  minimize the prediction error will be $c_k= -a_k$ of Eq(5).


how can I make a simple lowpass filter with this constraints?

check out the music-dsp source code archive here: http://www.musicdsp.org/showmany.php

I recommend a "Chamberlin" State-Variable Filter (which is 12dB/oct) which gives you other flavors (highpass, bandpass, etc) for free. If you cascade two of them you can get 24dB/oct, etc.


Basic difference between Multi-scale and Multi-resolution images

In scale space the square of width of blurring kernel considered as the scale (variance of Gaussian and not its Standard deviation because if you sequentially blur an image with 2 Gaussian kernel variance of effective blurring kernel will be sum of variance of each kernel) but in multi-resolution analysis the resolution defined by the scaling parameter.


Integral of an upsampled signal, without actually resampling it

If you upsample a causal discrete-time signal $x[n]$ by an integer factor $L$ and you use an interpolation filter with a causal impulse response $h[n]$, the upsampled and interpolated signal is given by

$$y[n]=\sum_{k=0}^{\infty}x[k]h[n-kL]\tag{1}$$

(Note that the requirement of causality can be easily removed; I just use it to have the lower summation indices start with the value $0$.)

If you want to compute the sum of $y[n]$, you would need to compute

$$s_y=\sum_{n=0}^{\infty}y[n]=\sum_{n=0}^{\infty}\sum_{k=0}^{\infty}x[k]h[n-kL]\tag{2}$$

However, you can interchange the sums in (2), resulting in

$$s_y=\sum_{k=0}^{\infty}x[k]\sum_{n=0}^{\infty}h[n-kL]\tag{3}$$

The value of the inner sum in (3) is independent of the index $k$, so we can write

$$s_y=s_h\cdot \sum_{k=0}^{\infty}x[k]=s_h\cdot s_x\tag{4}$$

with

$$s_h=\sum_{n=0}^{\infty}h[n]$$

So the sum of the upsampled and interpolated signal is simply the sum of the original signal times a factor which is given by the sum over the impulse response values of the interpolation filter. Note that for the ideal sinc filter, this factor equals the upsampling factor $L$. This is also at least approximately true for all practical interpolation filters, because the sum over the impulse response equals the DC gain of the interpolation filter, which should ideally equal $L$.

So if you want to keep things simple (and why wouldn't you?), you could use

$$s_y\approx L\cdot s_x$$

which is a very good approximation (or even exactly true) for any reasonable interpolation filter.

In the above derivations I've assumed that all the (infinite) sums exist, which in practice is always the case if we deal with finite length signals and with stable interpolation filters.