The plot has the power already measured for each sample. If the raw data samples (and not the power results shown here, then the average power can be direclty measured (given by the variance of the waveform). We may think that simply averaging the power measurements in the plot would give us the average signal power--- This would be the average of the squares which is not the same as a true power measurement given by the square of the average. It will be close but such an approach is highly dependent on the statistics of the waveform. For example, if the waveform was a complex Gaussian signal, the average of the squares would under-estimate the true power by -1.05 dB. This is detailed nicely in the app note by Agilent (formerly HP) # 1303 available at this link http://hpmemoryproject.org/an/pdf/an_1303.pdf .
Typically SNR would be determined from the noise power within the occupied bandwidth of the signal while a white noise contribution is spread evenly over the sampling bandwidth. If the signal is oversampled, then you would want to increase the white noise source accordingly to maintain a target spectral density within the bandwidth of the signal.


What is the result of fourier transform of an image

A 2D FFT consists in comparing 2D data with a discrete set of evenly spaced 2D complex exponentials $e_{u,v}$ (made of a cosine for the real part, and and a sine for the imaginary part). Behind the scene, this comparison is exact: a comparison of a 2D data with enough 2D complex exponentials bears the same information as the original data, albeit easier to interpret.

For each complex exponential (as numerous as pixels in the image), this comparison (or correlation) yields one complex coefficient $c_{u,v}$. The modulus of $c_{u,v}$ tells you about the proportion of amplitudes related to complex exponential $e_{u,v}$; the phase provides information on how the complex exponential contribution is shifted in the image. All in all, this provides you with a novel representation of the 2D data, as a sum of shifted exponentials with different amplitudes.

So, the complex values are "complex amplitudes", corresponding to frequency values indexed by the indices $_{u,v}$.


How to prove: high-slope/more change is high freq part of signal

The Fourier transform of the derivative of a signal $x(t)$ is given by

$$\dot{x}(t)\Longleftrightarrow j\omega X(\omega)\tag{1}$$

where $\dot{x}(t)$ denotes the derivative of $x(t)$, and $X(\omega)$ is the Fourier transform of $x(t)$. From (1) we have

$$\dot{x}(t)=\frac{j}{2\pi}\int_{-\infty}^{\infty}\omega X(\omega)e^{j\omega t}\,d\omega\tag{2}$$

If we now assume that $x(t)$ is band-limited with cut-off frequency $\omega_c$ (in radians), then (2) becomes

$$\dot{x}(t)=\frac{j}{2\pi}\int_{-\omega_c}^{\omega_c}\omega X(\omega)e^{j\omega t}\,d\omega\tag{3}$$

From (3) we can determine an upper bound for $|\dot{x}(t)|$:

$$|\dot{x}(t)|\le \frac{1}{2\pi}\int_{-\omega_c}^{\omega_c}|\omega| |X(\omega)|\,d\omega
\le\frac{\omega_c}{2\pi}\int_{-\omega_c}^{\omega_c}|X(\omega)|\,d\omega\tag{4}$$

The upper bound (4) on the magnitude of the derivative of $x(t)$ increases with increasing $\omega_c$, so you can expect faster changes of $x(t)$ with increasing values of the cut-off frequency $\omega_c$.


Kalman filter for position and velocity: introducing speed estimates

Is this a valid way of introducing speed estimates into the process?

If you choose your state appropriately, then the speed estimates come "for free". See the derivation of the signal model below (for the simple 1-D case we've been looking at).

Signal Model, Take 2

So, we really need to agree on a signal model before we can move this forward.  From your edit, it looks like your model of the position, $x_k$, is:

$$
\begin{array}
xx_{k+1} &=& x_{k} + \dot{x}_{k} \Delta t + \frac{1}{2} a (\Delta t)^2\\
\dot{x}_{k+1} &=& \dot{x}_{k} + a \Delta t
\end{array}
$$

If our state is as before:
$$
\begin{align*}
\mathbf{x}_{k}  & =\left(\begin{array}[c]{c}x_{k}\\
\dot{x}_{k}\end{array}
\right)
\end{align*}
$$
then the state update equation is just:
$$
\mathbf{x}_{k+1} =
\left(\begin{array}[c]{c} 1\ \ \Delta t\\
0\ \ 1\end{array}
\right)
\mathbf{x}_{k} +
\left(\begin{array}[c]{c} \frac{(\Delta t)^2}{2} \\
\Delta t \end{array}
\right)
a_k
$$
where now our $a_k$ is the normally distributed acceleration.

That gives different $\mathbf{G}$ matrix from the previous version, but the $\mathbf{F}$ and $\mathbf{H}$ matrices should be the same.

If I implement this in scilab (sorry, no access to matlab), it looks like:

// Signal Model
DeltaT = 0.1;
F = [1 DeltaT; 0 1];
G = [DeltaT^2/2; DeltaT];
H = [1 0];

x0 = [0;0];
sigma_a = 0.1;

Q = sigma_a^2;
R = 0.1;

N = 1000;

a = rand(1,N,"normal")*sigma_a;

x_truth(:,1) = x0;
for t=1:N,
    x_truth(:,t+1) = F*x_truth(:,t) + G*a(t);
    y(t) = H*x_truth(:,t) + rand(1,1,"normal")*sqrt(R);
end

Then, I can apply the Kalman filter equations to this $y$ (the noisy measurements).

// Kalman Filter
p0 = 100*eye(2,2);

xx(:,1) = x0;
pp = p0;
pp_norm(1) = norm(pp);
for t=1:N,
    [x1,p1,x,p] = kalm(y(t),xx(:,t),pp,F,G,H,Q,R);
    xx(:,t+1) = x1;
    pp = p1;
    pp_norm(t+1) = norm(pp);
end

So we have our noisy measurements $y$, and we've applied the Kalman filter to them and used the same signal model to generate $y$ as we do to apply the Kalman filter (a pretty big assumption, sometimes!).

Then following plots show the result.

Plot 1: $y$ and $x_k$ versus time.

Plot 2: A zoomed view of the first few samples:

Plot 3: Something you never get in real life, the true position vs the state estimate of the position.

Plot 4: Something you also never get in real life, the true velocity vs the state estimate of the velocity.

Plot 5: The norm of the state covariance matrix (something you should always monitor in real life!).  Note that it very quickly goes from its initial very large value to something very small, so I've only shown the first few samples.

Plot 6: Plots of the error between the true position and velocity and their estimates.

If you study the case where the position measurements are exact, then you find that the Kalman udpate equations produce exact results for BOTH position and speed.  Mathematically it's straightforward to see why.  Using the same notation as the wikipedia article, exact measurements mean that $\mathbf{z}_{k+1}=x_{k+1}$. If you assume that the initial position and speed are known so that $\mathbf{P}_k=0$, then $\mathbf{P}_{k+1}^{-}=\mathbf{Q}$ and the Kalman gain matrix $\mathbf{K}_{k+1}$ is given by

$$
\mathbf{K}_{k+1} = \left(\begin{array}[c]{c}1\\
2/dt\end{array}
\right)
$$

This means that the Kalman update procedure produces

$$
\begin{align*}
\mathbf{\hat{x}}_{k+1} & = \mathbf{F}_{k+1}\mathbf{x}_k + \mathbf{K}_{k+1}\left(\mathbf{z}_{k+1} - \mathbf{H}_{k+1} \mathbf{F}_{k+1}\mathbf{x}_k\right)\\
& = \left(\begin{array}[c]{c}x_k + \dot{x}_k dt\\
\dot{x}_k\end{array}
\right) +
\left(\begin{array}[c]{c}1\\
2/dt\end{array}
\right)
\left(x_{k+1} - \left( x_k + \dot{x}_k dt\right)
\right)\\
& = \left(\begin{array}[c]{c}x_{k+1}\\
 2 \left(x_{k+1} - x_k \right) /dt - \dot{x}_k\end{array}
\right)
\end{align*}
$$

As you can see, the value for the speed is given by exactly the formula you were proposing to use for the speed estimate.  So although you couldn't see any kind of calculation $(x_k-x_{k-1})/dt$ for speed, in fact it is hidden in there after all.


How to Count the number of People in a room using CCTV Camera feed?

This is clearly not an easy task.

The problem is, if you want a more-or-less accurate count, then you need to turn to advanced algorithms (and maybe use 2 cameras, or a kinect).

If you can't afford to take this path, then you need to try simpler options. Personally, I would try the following:

detecting skin pixels, segmenting the image with respect to that, then count the blobs (complexity medium);
detect faces, then count them. Obviously, people have to be facing the camera, but you can find any implementations of face detection with different libraries (complexity low);
use background subtraction (after some careful tuning to take into account natural illumination changes) and guess a count based on the number of pixels that belong to the foreground (complexity low, but expect poor results);
apply some existing pedestrian detection algorithm (complexity medium to high).

You don't have to implement all that by yourself. In your case, checking what OpenCV can do (or if your algorithm of choice is implemented in OpenCV) seems a good choice because:

Processing and openFrameworks have wrappers to OpenCV;
OpenCV usually has a good runtime performance;
many computer vision algorithms that will interest you are developed with OpenCV (face detection, background subtraction, probably pedestrian detection).


Statistics to characterize time-domain envelope of acoustic signal

An ADSR has exactly five parameters, and you need these five to characterize it - every set of parameters yields a different envelope.
And frankly, the ADSR envelope is pretty easy to visualize, and to understand:
So, calculate the envelope of your signal, find the

peak amplitude,
the time it takes to reach that peak
the plateau level
the time it takes to reach that plateau from the peak, and
the time it takes for the envelop to decay after leaving the plateau.

These are your five degrees of freedom, and they can be quite directly read from a slightly averaged envelope.


How to estimate the mean of a clipped signal?

Is this crazy? Good? Done before and has a name?

I wouldn't say that it is "Good" because clipping is not noise that can be modeled as an additive component with some distribution.

  fit a curve over the top 3 or 5 highest bins, and the position of the
  peak is the signal median

The position of that peak would be the modal mean, not the median

The effect of clipping on the signal's distribution would be to a) truncate the real distribution (whatever that is) at the maximum (which depends on the sensor & measurement (or "environment" as you put it) and b) distort the distribution towards this maximum depending on how much clipping was present in the signal.

b) means that whatever you do, your median / mean will be biased by the amount of clipping.

Ideally, you would have to change the sensor to one that can respond to higher accelerations.

If you absolutely must use the current sensor and deal with the clipping then prior to any further processing of its output you could try removing the clipped values (with a simple thresholding operation) and fitting a curve to the data. If you were to fit an spline for example, then the curve would depart at the beginning of the clipped part with the slope of the usable signal prior the clipping and arrive at the end of the clipped part with a slope of the usable signal after the clipping. In this way, all clipped parts will be converted to some "soft" peak (because of the splines "smoothing" effect).

Of course, this is a compromise. It could well be that the clipping is so much that very little "usable" signal remains on either side to fit a "good" spline. Furthermore, depending on your sampling frequency ($Fs$), the signal might be full of impulse-like data which would lead to "very tall" splines as the algorithm would try to fit the sharp slopes. More information would be needed to take a decision on that but in general a higher $Fs$ would help.


How to compute the channel gain from the path loss model of a wireless channel

To answer the literal questions you asked:

yes
see above

To Frame Challenge your question and give you what you were probably looking for:

The path loss equation you cited:
$$
128+37\log_{10}(d) ,
$$
is equivalent to:
$$
\frac{P_r}{P_t} = \frac{10^{-12.8}}{d^{3.7}},
$$
where $P_r$ is the received power in Watts, and $P_t$ is the transmit power in Watts.  From this equation it can be seen that the path loss exponent is 3.7.

Under log-normal shadowing, the random shadowing component is modeled as a zero-mean Gaussian random variable (in dB) added to the path loss.


Understanding and Correcting phase shift in General Linear Phase FIR filters

A generalized linear phase filter with $\beta=\pi/2$ or $\beta=3\pi/2$ introduces (apart from a constant delay) a frequency independent phase shift of $\pm 90$ degrees. Two standard filters that do exactly that are a differentiator with frequency response

$$H(\omega)=j\omega\tag{1}$$

and a Hilbert transformer with frequency response

$$H(\omega)=-j\,\text{sign}(\omega)\tag{2}$$

Note that both frequency responses are purely imaginary and conjugate symmetric, i.e. $H(\omega)=H^*(-\omega)$, as required for real-valued systems.

In the time domain, the filters do the following:

$$y(t)=\frac{dx(t)}{dt}\qquad\text{(differentiator)}\\
y(t)=\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{x(\tau)}{t-\tau}d\tau\qquad\text{(Hilbert transformer)}$$

So in general these systems introduce severe distortion. It is relatively easy to imagine the effect of differentiation. To get a feeling for what the Hilbert transform does, take a look at this table showing some transform pairs.

These two systems can be inverted. The differentiator can be inverted by an integrator, and the Hilbert transformer can be inverted by another Hilbert transformer and by sign reversal.

Note that these filters are just two common examples of systems with imaginary frequency responses, but they are representative in the sense that they introduce severe distortions in the time domain waveform, just like all other systems with imaginary frequency responses with conjugate symmetry.


Accelerometer Inference from power spectrum

What could the cause of this trend be? System bias or something else?

Most mechanical systems are low-pass, i.e. they let through less energy at higher frequencies. That would be perfectly plausible, but I don't know what you're observing.
Of course, your accelerometer is inherently frequency-selective, too, and you might need to calibrate it, together with your recording setup.
Lastly, we don't know exactly how you post-process.

Best way of handling this - Matlabs detrend tool does not have much of an effect. What would be the best way to correct this?

"Correcting  this" means "I saw something that is clearly visible in my measurement, but I don't know where it comes from and whether it's actually there in the observed phenomenon. Let me get rid of the evidence".
That's called a lie.
So, before changing anything about the symptom (plot), make sure you're not masking a cause. This is less about the signal processing than about physically understanding your measurement setup. It can, and probably will, involve you measuring known acceleration amplitudes at different frequencies, and looking at whether the results are as expected, or whether you need multiply with a different factor for different frequencies.


Removing marker pen artifacts from scanned microscope slide images

You could also try local contrast enhancement or similar. Fiji has a function called Enhance Local Contrast (CLAHE) which does this. Looks like high pass filtering works reasonably well for you though.

As for your second question. The Gaussian Blur wikipedia you mention refers to Low pass filter. For the low pass filter each pixel is average over the surrounding pixels, weighted by a Gaussian centred on the pixel. In your case you are doing a high pass filter. This can can be given by the original minus the low pass filter, which is probably the easiest way of thinking about it (IMO).


Counter-example to show that this system is non-linear?

OK. Let's try homogeneity again:
$$
y[n] = x[n]\cdot x[n-1]
$$
for input $x[n]$.

For input $a\cdot x[n]$ we get
$$
y'[n] = a\cdot x[n] \cdot a x[n-1] = a^2 \cdot x[n]\cdot x[n-1] \not = ay[n]
$$
so homogeneity doesn't apply.

Let's try that with an example, as per MBaz's suggestion.

Let's try:
$$
x[n] = 2
$$
so that
$$
y[n] = 4
$$
Then let's try twice that:
$$
x'[n] = 2\cdot x[n] = 4
$$
So that
$$
y'[n] = 16 \not= 2\cdot y[n]
$$


High pass filter of 0.01 Hz for EEG slow drifts

You're asking for a time machine: in order to avoid the DC bump at $t=0$, you need to know the state of the filter as if it had been running before the recording ever began. That's not going to happen no matter what kind of filter you use.

Here's a trick which may come in handy for you. Let's say you have an $N$-tap FIR filter (or an IIR filter whose impulse response decays adequately to 0 after $N$ samples). Take the first $N$ samples of your signal, reverse them, and prepend that to the start of your signal. In effect, we're defining a new signal $g(t)=g(-t)=f(t)$ for the purposes of filtering. When you filter that, the filter state will have already been "primed" with the initial signal: delete the first $N$ samples of the output, and the DC bump should have hopefully gone away.

Needless to say, there are numerous gotchas with that scheme — for starters, it's harder to do for realtime processing, and if the signal derivative is high at $t=0$, you are still going to have trust issues with the first few samples. But it should get rid of the bump.


What are best practices to compute an audio spectrogram?

I think, perhaps, that the statement should be reversed: that a lot of implementations "just" use the square magnitude. This is because the DFT of a real-valued signal is typically complex, so that the square magnitude is just the product of complex conjugates, e.g.:

  X*(f) · X(f)

This requires less computation than the magnitude, which would follow the same algorithm, but with a subsequent square-root.

The square magnitude of a DFT is formally known as the Power Spectrum.


this is what we sometimes call "Vector Synthesis" and other times call "Wavetable Synthesis".  linear crossfading between equal-sized circular lookup tables usually suffices.

but make sure the waveforms are aligned.  you can circularly delay (or "rotate") one lookup table relative to the other.  choose a relative rotation so that the cross-correlation between the two waveforms is maximum.  if you somehow crossfade between a sine wave and a square wave where the fundamental of the square wave is outa phase with the sine, you'll get a glitch in the crossfade.


Guidance/Brainstorming for a mapping/classification problem

You need to be certain that your feature vector is capturing the properties you are looking for. Using the mean and median for dimension reduction (i.e. extracting the features) is not a good choice unless your problem is really simple. Better choices are: PCA, ICA, and SVM's. Capturing the features in the frequency domain can also be a great tool if your signal is sparse in the FFT or DCT domain.


Some questions about the intuition of the DWT

First, one should be cautious about processing short signals like these. I am unsure about the length of 3 for a3 and d3. Now, the many questions:

I would not call them "frequency plot", but stacked subband plots. They are just illustrations of the behavior of the wavelet coefficients. For each subband, different scalings may be used. Another rendering is stacking the intensity of coefficients in a image. For the interpretation of the points: the coefficients in each subband at located in time and scale, so the interpretation in $[\pi/2T, \pi/T]$, $[\pi/4T, \pi/2T]$ is only approximative (because of aliasing, subsampling, etc.).
No, one is concatenating chunks of coefficients in different subbands and at different rates
Yes, subbands of a given WPT can generally be interpreted in dyadic portions of the time-scale axis
Haar wavelets (2-tap filters) have the rare property (among standard wavelet filters) they don't overlap at a given scale. For longer filters with 2-band wavelets, those filters overlap, and yield troubles at the extremities of a finite-length signal. There is a lot of works on non-expensive expansions, relaying on the symmetry/antisymmetry of the wavelet.


This is how I would approach the problem.

System $A$ cannot be LTI because its output has frequencies not present in its input. The input spectrum is a brickwall from $f=-1/2$ to $f=1/2$, but the output spectrum is a brickwall from $f=-1$ to $f=1$.

It's impossible to arrive at a conclusion for system B. It could be an ideal LTI low-pass filter, but it could also be something else (for example, it could be a signal generator).

Note that we cannot say that A and B are inverses of each other. That statement requires $x(t)=B[A[x(t)]]$ to be true for all $x(t)$. Here, we only know that is true for one specific input.


Is $\arctan(t)$ an energy or power signal?

The typical inverse tangent function maps the input range of $ t \in (-\infty,\infty)$ into an output range of $(-\pi/2, \pi/2)$ as in the figure below:

Based on this, its values are bounded for all $t$. Yet, since the intergal of its square is unbounded, then it cannot be an energy signal;i.e.,

$$ \int_{-\infty}^{\infty} |\tan^{-1}(t)|^2 dt ~~~~~\to \infty $$

Moreover the average power of this signal is nonzero and finite, hence this is a power signal;i.e.,

$$ 0 < \lim_{T \to \infty} \frac{1}{T}\int_{-T/2}^{T/2} |\tan^{-1}(t)|^2 dt ~~~~~ < \infty $$

To evaluate the last integral, use the following MATLAB figure to see that

:

the square of the inverse tangent is always less than $(\pi/2)^2$ for all $t$. Hence to find an upper bound on the integral we can use $|\tan^{-1}(t)|^2 < (\pi/2)^2 = K$. Then the integral becomes

$$ 0 < \lim_{T \to \infty} \frac{1}{T}\int_{-T/2}^{T/2} |\tan^{-1}(t)|^2 dt < \lim_{T \to \infty} \frac{1}{T}\int_{-T/2}^{T/2} K dt = \lim_{T \to \infty} \frac{1}{T} K \cdot T = K ~~~~~ < \infty $$


From a 44kHz audio signal, how could I extract 32 evenly spaced frequency bands?

Depends on what you want to do with it. Some possibilities:

Short Time Fourier Transform with a 64 FFT length will give you 33 bins as a function of time. These are actually 31 "normal" bins (with amplitude and phase) and 2 "half" bins that only have amplitude (DC and Nyquist) and half the bandwidth. If your input vector is complex, you can use an FFT length of 32 and you get 32 "normal" bins.
Do an FFT over the whole sequence, create any bins you want in the frequency domain and simply integrate the energy in each bin. However this will discard all phase information.
If you really need 32 bit equally spaced complex frequency bins for
a real valued input signal you can do "half bin rotation" of the
frequency sampling grid. You basically multiply your sequence with a
exp(-jw/N * (n/2)) before the FFT grid. YOur first bin would then go from 0Hz to 1378Hz instead from -689Hz to 689 Hz.

It really depends on how you want to do your frequency binning. Do you care about the location of center frequencies or the band edges? Do you need sharp edges or overlaps? Do you need amplitude and phase? Is the input real or complex? All this factors into the correct binning procedure.

Let's look at the example of calculating the A-weighted sound pressure level from a time domain microphone signal. Your input is real valued, i.e. you can ignore the negative frequencies. You care about center frequencies but you need a log spaced grid. You want some overlap between the bands so that frequencies that are close to a band edge have some contribution in both bands. However the overlap should conserve energy, so the overall energy is frequency independent. We don't care about phase, so we can simply integrate energy inside a band to get a single number.

As this example shows, frequency binning is application dependent on not as trivial as one might initially think.


Audio Signal Separation - Identifying Interest points

(a follow-up to my suggestion on the previous question), you can use the spectrogram and ICA to help:

A similar shorter sound file:

import wave, struct, numpy as np, matplotlib.mlab as mlab, pylab as pl
def wavToArr(wavefile):
    w = wave.open(wavefile,"rb")
    p = w.getparams()
    s = w.readframes(p[3])
    w.close()
    sd = np.fromstring(s, np.int16)
    return sd,p

def wavToSpec(wavefile,log=False,norm=False):
    wavArr,wavParams = wavToArr(wavefile)
    print wavParams
    return  mlab.specgram(wavArr, NFFT=256,Fs=wavParams[2],detrend=detrend_mean,window=window_hanning,noverlap=128,sides='onesided',scale_by_freq=True)

wavArr,wavParams = wavToArr("bat_speech.wav")
hf = pl.figure(); ax=hf.add_subplot(1,1,1)
ax.plot(wavArr)

Now take a look at the spectrogram:

Pxx, freqs, bins = wavToSpec("bat_speech.wav")
Pxx += 0.0001
freqs += (len(wavArr) / wavParams[2]) / 2.
hf=pl.figure(figsize=(12,12));
ax = hf.add_subplot(2,1,1);
#plot spectrogram as decibals
hm = ax.imshow(10*np.log10(Pxx),interpolation='nearest',origin='lower',aspect='auto')
hf.colorbar(hm)
ylcnt = len(ax.get_yticklabels())
ycnt = len(freqs)
ylstep = int(ycnt / ylcnt)
ax.set_yticklabels([ int(freqs[f]) for f in xrange(0,ycnt,ylstep) ])

We can clip this at 8000Hz or so it looks like, or don't bother cleaning it up.

Now you have frequencies which can proxy for multiple sources in BSS, so you can play with PCA, ICA, normalization, etc.  For example, see if you have some components you can isolate:

from sklearn.decomposition import PCA, FastICA
ncomps = 7
# reduce dimensionality with PCA
pca = PCA(n_components=ncomps)
y = Pxx.copy().T
pc = pca.fit(y).transform(y)
# run ICA
ica = FastICA(n_components=ncomps,random_state=42)
z = ica.fit(pc).transform(pc).T
hf = pl.figure()
for p in xrange(ncomps):
    ax = hf.add_subplot(ncomps,1,p+1)
    ax.plot(z[p])

or see if the spectrogram is enough to let you do your segmentation:

hf = pl.figure()
ax = hf.add_subplot(1,1,1)
ax.plot(np.sum(Pxx,axis=0))

EDIT:  Just realized you said you didn't want to use the frequency domain, but it may help you isolate your phonemes which you can extract.  Anyway, running on your SeroWeb.wav I get this:


frequency folding due to aliasing in FFT

The Fourier transform $X(e^{j\Omega})$ of a sequence $x[n]$ is by definition periodic of period $2\pi$ in $\Omega$. The $N$-point discrete Fourier transform (DFT) of $x[n]$ is simply defined as a sampled version of the Fourier transform $X(e^{j\Omega})$, that is
$$X[k] = X(e^{j\frac{2k\pi}{N}}).$$
The periodicity of $2\pi$ in $\Omega$ for the Fourier transform becomes a periodicity of $N$ in $k$ for the DFT
$$X[k+N] = X[k].$$
The negative frequencies showing up to the right of the positive frequencies is simply a consequence of this periodicity. Or if you prefer, you get $X[k]$ for $k \in [0, N-1]$ instead of $k \in [-\frac{N-1}{2}, \frac{N-1}{2}]$. If you want the vector to contain first negative frequencies, then the DC and finally the positive frequencies instead, you can for example use MATLAB fftshift.

Those negative frequencies are not shown for the "real" DFT because if $x[n]$ is real, then $X[k] = X^*[-k]$ (conjugate symmetry). Hence, you only need $X[k]$ for $k \in [0, N/2]$.


I can't comment because I'm a new user... have a look to these questions:

Estimating the time offset (delay) between two audio signals in real-time

Generalized Cross Correlation

Time delay estimation of oscilloscope signals using cross correlation


Understanding BFSK Modulation (MATLAB/Simulink)

So, ignoring noise, the input into the receiver multiplier coming from the channel is either $\cos(2\pi f_0t)$ or $\cos(2\pi f_1t)$. The other multiplier input is $\cos(2\pi f_0t)-\cos(2\pi f_1t)$. All signals are defined from 0 to $T$. Let's assume the transmitted signal has frequency $f_0$. Then you have an integrator:

\begin{align}
&\int_0^T \cos(2\pi f_0t) \left[\cos(2\pi f_0t)-\cos(2\pi f_1t)\right]dt \\
&=\int_0^T \cos(2\pi f_0t)\cos(2\pi f_0t) dt - \int_0^T \cos(2\pi f_0t)\cos(2\pi f_1t)dt.
\end{align}

The frequencies $f_0$ and $f_1$ are chosen so that the cosines are orthogonal for $0<t<T$; that means that
$$\int_0^T \cos(2\pi f_0t)\cos(2\pi f_1t)dt =0.$$
The integrate and dump block outputs the result of the definite integral at $t=T$. You can see that the result is positive if frequency $f_0$ was transmitted and negative if $f_1$ was transmitted.

The sign block is a quantizer, and the lookup table converts the numbers back to bits 0 and 1.


Moving average algorithm that conserves integral

Why do you want to preserve the the sum? Is it to make the filter with the same gain everywhere?

If you want to keep the sum the same you need to pad the beginning and the end of the signal as:
$$
.... ,y_2,y_1,y_0 ---- y_0,y_1, y_2,...., y_n ---- y_n, y_{n-1},y_{n-2}...
$$
Example:
suppose that you have 6 samples

$$
y_0, y_1, y_2, y_3,y_4, y_5
$$

and you want to to do a moving average with length 5:

Special treatment at the beginning:

$$
\hat y_0   = (y_1  +  y_0 +  y_0 +  y_1 + y_2)/5,
$$
$$
\hat y_1   = (y_0  +  y_0 +  y_1 +  y_2 + y_3)/5
$$

Normal filtering:

$$
\hat y_2   = (y_0  +  y_1 +   y_2 + y_3  + y_4)/5,
$$
$$
\hat y_3   = (y_1  +  y_2 +   y_3 + y_4  + y_5)/5,
$$

Special treatment at the end:

$$
\hat y_4   = (y_2  +  y_3 +   y_4 + y_5  + y_5)/5,
$$
$$
\hat y_5   = (y_3  +  y_4 +   y_5 +  y_5 + y_4)/5,
$$


Linear response function for a system with derivative: $U=L \frac{d I}{dt}$, expressing $U=f(I)$

Let's start with the simple system

$$y(t)=x(t)\tag{1}$$

where $y(t)$ denotes the output and $x(t)$ denotes the input. It is linear and time-invariant, and, consequently, it must be possible to represent it by the convolution integral

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

where $h(t)$ is the system's impulse response. In this case, it turns out that $h(t)=\delta(t)$, i.e., the impulse response is a Dirac delta impulse.

Now, for a differentiator

$$y(t)=x'(t)\tag{3}$$

we "simply" (conceptually, at least) have $h(t)=\delta'(t)$, where $\delta'(t)$ is the (generalized) derivative of $\delta(t)$. So, using the concept of a generalized derivative of a distribution, we can write the input-output relationship of an ideal differentiator as the following convolution integral:

$$y(t)=\int_{-\infty}^{\infty}x(\tau)\delta'(t-\tau)d\tau\tag{4}$$


How to implement a filter associated to a specific wavelet


What you have is a very sharp rectangular boxcar shaped function.  It is approximately centered at the mid point of your data set.  Your magnitude spectrum looks correct, in that it is a really low fundamentally frequency with a ton of harmonics.  Your FFT without ‘absolute values’ is presumably plotting just either the real or imaginary values.  The oscillation in that plot occurs because of where that pulse is centered in the data set.  It’s unlikely that plot has any useful information.  I have no idea what you are trying to find from this endevour, but if you were interested in the period of the pulse you could zero of the end of the data set, providing more granular frequency data that you could use to measure the fundamental.  That or you could just measure it directly, which is probably a lot easier.


How to extract noise when the noise on three data sets are near identical?

Given the three measurements $y1$, $y2$ and $y3$ coming from three unknown signals $x1$, $x2$ and $x3$, which we could assume additive to some unknown common noise $x4$, we can write:
$$
y1=x1+x4\\
y2=x2+x4\\
y3=x3+x4
$$
or equivalently:
$$
y=Ax\\
A=[1,0,0,1;0,1,0,1;0,0,1,1]
$$
Again, assuming a enough noisy signal, we can have an standard Recursive Least Squares algorithm for solving the standard LS projection recursively:
$$
\hat x=(A^TA)^{-1}A^Ty
$$
This way, and if and only if the additive noise and the noisy signal suficciencies are met, you will solve for your three signals, and your coupled noise.

Because of the former, you even can express this solution using the Kalman Filter


If an autocorrelation function of a process is a rectangular function, then is the process deterministic?

An autocorrelation function like that is not possible for ordinary signals. The autocorrelation function is the Fourier transform of the power spectral density, which is a strictly positive quantity. However, to get a rectangular function as Fourier transform you need to have negative values. That is in conflict with the strictly positive power spectrum.


Determining Loss of Information by Taking Average (Mean of Signal)

Consider observing $y(n) = w(n)$ where $w(n)$ are IID gausian distributed with mean 0 and variance $\sigma^2$. The entropy of each other observations $y(n)$ is thus: $$\log_22\pi e\sigma^2$$.

Now if you average $y(n)$ the resultant variable is also gaussian distributed with mean 0 but variance $\sigma^2/N$. Therefore the entropy is $$\log_22\pi e\sigma^2/N$$ so the entropy difference is $$\log_2N$$ and this would be the entropy reduction.

The amount of information lost will depend on the distribution of the random variable. However, there are some interesting problems in signal processing where there would essentially be no loss.

For ex: consider a set of observations taken as $$y(n) = x + w(n)$$ where $w(n)$ is normally distributed IID random varibles and suppose we take $N$ observations of this kind.

Then the mean of $y(n)$ essentially represents what we refer to as a sufficient statistic for $x$ and infact represents the minimum variance unbiased estimator, i.e. the most efficient estimator. So in this case the sample mean encompasses all the necessary information and no information is lost.

However, if the noise $w(n)$ were  not gaussian we couldn't claim this. So as you see the amount of information lost by sample mean depends on the distribution of the underlying variable.

Another example would be a trivial one of a deterministic random varaible. Here also no information would be lost by the sample mean.

Another way to look at it would be if the likelihood of posterior distribution (parameterized by the sample mean), after observing the sample mean renders the posterior useless, meaning encompasses no additional information about the random variable then there is no information loss, otherwise we would incur information loss.

Hope that helps


FFT of color images incorporated into an Object Recognition method?

Okay I found some ideas about using Fourier Transforms in Object Recognition. What I am unable to understand is however, its downsides. Why was the method discarded.

For those who are interested in reading about how Fourier Transforms were applied in Object Recognition , look for Fourier Mellin Transforms.

The Fourier-Mellin transform is a useful mathematical tool for image recognition because its resulting spectrum is invariant in rotation, translation and scale. The Fourier Transform itself (FT) is translation invariant and its conversion to log-polar coordinates converts the scale and rotation differences to vertical and horizontal offsets that can be measured. A second FFT, called the Mellin transform (MT) gives a transform-space image that is invariant to translation, rotation and scale.

The steps:
The Fourier Transform

The Discrete Fourier Transform (DFT) is given by the following expression:

This is often computed using the Fast Fourier Transform algorithm to speed up the time needed for the calculation.

Cartesian to Log-Polar conversion
The FFT is projected onto the log-polar plane by the coordinate transform shown below.

For the conversion from Cartesian coordinates to Log-Polar coordinates the following equation is true:

$$r = \sqrt{x^2+y^2}$$

The origin $(m_o, n_o)$ should be at the center of the image matrix to ensure the maximum number of pixels is included. If the image consists of a square $N \times N$ matrix then the coordinates of the origin are:

The maximum sampling radius for the conversion can now be calculated as:

If an inscribed circle is chosen as the conversion boundary, some pixels that lie outside the circle will be ignored. If a circumscribed circle is chosen, all pixels will be taken in account, but some invalid pixels will also be included (pixels falling inside the circle but outside of the image matrix).

Since the pixels in Cartesian coordinates cannot be mapped one-to-one onto pixels in the Log-Polar coordinate space, an average of the surrounding pixels needs to be calculated. The standard methods to do this includes nearest neighbor, bilinear and bicubic resampling.

The relationship between the polar coordinates ($\rho$, $\theta$) used to sample the input image and the polar coordinates of the log-polar image (r , $\theta$) can be described by:

$$(\rho,\theta)    =  (e^r, \theta) $$

To map the input image pixels imagein(xi,yi) onto the output image pixels imageout(rm, $\theta$m), the co-ordinates xi, yi are computed using:

where $$(\rho,\theta)    =  (e^r, \theta) $$

The next step is to get a transform-space image that is a rotation and scale invariant representation of the original image. This is done using the Mellin transform which can be expressed as:

Converting to polar coordinates, we have:

Now using $$(\rho,\theta)    =  (e^r, \theta) $$

we get:

The Fourier Mellin transform has been used for Image Recognition in R. Moller, H. Salguero, E. Salguero: “Image recognition using the Fourier-Mellin Transform”, LIPSE-SEPI-ESIME-IPN, Mexico


I have done this before and it works as long as you tune your parameters correctly. There should be a gain factor on the integrator which is < 1.0.

If your application is audio, one thing to watch for is audible pitch changes during periods when the ratio is changing. You can add a proportional term to your loop filter to tune the transient response if this is an issue.


There may be some confusion about the definition of SNR.

In some applications you may be interested in reducing the ratio of the signal power spectral density to noise power spectral density.  "SNR" in this case refers to the ratio of the signal to noise power when you measure both over the same bandwidth.  In these cases, as Juancho pointed out, direct-sequence spread spectrum is a classical technique.  It deliberately spreads the signal power over a wider bandwidth.

The SNR used in the Shannon-Hartley Theorem refers to the SNR per channel usage.  To oversimplify the matter, it refers to the SNR per symbol after you demodulate the signal.  Spread spectrum does not decrease the SNR needed to achieve a given data rate; it just spreads the signal power in the SNR calculation over a wider bandwidth.  It is true that you can communicate with SNR lower than 0 dB: one typically does so in a practical environment using simple modulation (e.g., BPSK, QPSK) and a low code rate with a strong-performing code.  Whether one spreads the signal or not is a separate decision.


I think that what you are looking for is a maximum intensity projection (MIP).
The general approach has two parts:

Define the view
Obtain the projection

To define the view, you need to either rotate the data in-place or rotate the view plane around the data and then use the coefficients of the plane to read through your data matrix. In either case, you are going to need interpolation.
The most straightforward way to do this is to rotate the data in-place.
If you are interested in obtaining just the three principal views, you don't even need to rotate the data.
Let's assume that you have a matrix like $(y,x,z)$. That is, each axial image's data is in $(y,x)$ and each slice's axial image can be accessed via $z$.
Since these are the axial images, looking at an individual $(y,x)$ image is like looking through our volume from the top. So, we are looking at $(y,x)$ down through $z$.
To obtain an "X-Ray-like" image of the volume as is, all you have to do is iterate through each $(y,x)$ pixel, get all the pixel values that "live" in that column (so, get all the pixel values at location $(y,x)$ from slice $0,1,2,3,..,N$) and select the maximum value from that column. Then move to the next pixel and so on.
This will give you the MIP in the axial direction. To obtain the MIP from another point of view, all you have to do is to rotate your volume to align it with the view and run the same iteration through the image again.
In the special case where you only need the three principal directions, you don't even need to rotate the data, you just change the way you are reading it: Where before we were looking at $(y,x)$ down through $z$ for the axial view, we now need to look at $(y,z)$ down through $x$ for the sagital view...and so on.
Here is a quick example in Python (but generalisable easily to other platforms too):
import numpy
from matplotlib import pyplot as plt

# We need some data, so let's build a phantom.
# This is a sphere that is harded in the middle and contains
# a rectangular block.

v = numpy.zeros((10,10,10)) # This is the size of our phantom

for m in range(0,10):
    for n in range(0,10):
        for k in range(0,10):
            # This is basically the a reversed distance plus a dense rectangular block running through it.
            v[m,n,k]= 20-numpy.sqrt((m-5)**2+(n-5)**2+(k-5)**2) + (10 if m>3 and n>3 and k>2 and m<6 and n<6 and k<8 else 0)

# To get the axial MIP:

z = numpy.zeros((10,10)) # Notice this is a single image, NOT a stack

for m in range(0,10):
    for n in range(0,10):
        z[m,n]=numpy.max(v[m,n,:])

plt.imshow(z);plt.show()

And this will give you something like:

This is a MIP of the axial view, so looking at our volume from the top.
To look at our volume "from the side" (the sagital view), we just change the way we read the volume:
# To get a sagital view MIP:

for m in range(0,10):
    for n in range(0,10):
        z[m,n]=numpy.max(v[m,:,n])

plt.imshow(z);plt.show()

And this will give you something like:

Notice here that the $x$ axis of this image is the $z$ axis of the volume.
To obtain any other view, simply 3D rotate the data in your v to align it to the view plane and then run the z iteration again for the MIP.
You can see here that if the pixel size and slice thickness are not equal, then the sagital view will appear distorted.
If you are looking for a library to achieve this, you can use something like vtk, although, be aware that it requires some level of prior experience to use it and more overhead to get an example up and running.
Hope this helps.


Obtaining impulse response by means of cross-correlation between input and output

we can calculate it!

define $y=h*t$, $r_{xy} -$
is the cross correlation

$r_{xy}(\tau)=\int_{-\infty}^ \infty x(t)y(\tau+t)dt = (y*x^-)(\tau) = (x^-*h*x)(\tau) =(r_{xx}*h)(\tau)\underbrace{=}_{r_{xx}=\delta}h(\tau)$

so actually the cross correlation is the impulse response


Active noise cancellation using kalman filter

I guess you use the LMS and RLS filters because they are adaptive, and you can do the same thing with the Kalman filter, if you augment the state-vector with the parameters of the model. You then get a non-linear system, and you can solve that using any number of non-linear Kalman-filter implementations, the most common being the extended Kalman filter (EKF), then second most common the unscented Kalman filter (UKF), and then there are a fair few more exotic variants, such as the second order Kalman filter. The convergence properties of any of these filters have not been proved, but they seem to work really well in practice. The UKF is typically regarded as better than the older EKF, since it does not require the inversion of the Jacobian, and it behaves as a higher order approximation -- 3rd order for Gaussian noise I seem to remember (compared to first order for the EKF).

The EKF is perhaps a bit more straight forward to get started with. You use the EKF for parameter adaptation by modeling unknown parameters as Wiener processes. Consider the linear system
\begin{align*}
 \dot{x} &= A(\theta) x + B(\theta) u + w \\
 y &= C(\theta) x + v
\end{align*}
with unknown parameters $\theta \in \mathcal{R}^{p}$. By augmenting the state vector to include the unknown parameters, $\chi^{\text{T}} = [x^{\text{T}}, \theta^{\text{T}}]$, we obtain the non-linear system $\dot{\chi} = f(\chi,u) + w$, $y = h(\chi) + v$, where
\begin{align*}
 f(\chi,u) &=
 \left[ \begin{array}{c}
  A(\theta) x + B(\theta) u \\
  0
 \end{array} \right] \\
 h(\chi) &= C(\theta) x .
\end{align*}
The Jacobians $F$ and $H$ for this system are found as
\begin{align*}
 F &=
 \left[ \begin{array}{cc}
  A(\theta) & \frac{\partial}{\partial \theta}[A(\theta) x + B(\theta) u] \\
  0 & 0
 \end{array} \right]_{\hat{x}, \hat{\theta}} , \\
 H &=
 \left[ \begin{array}{cc}
  C(\theta) & \frac{\partial}{\partial \theta}[C(\theta) x]
 \end{array} \right]_{\hat{x}, \hat{\theta}} .
\end{align*}

The general implementation of the EKF and UKF, as well as for parameter adaption, can be found in several publication. I personally like "Optimal State Estimation: Kalman, H$_{\infty}$, and Nonlinear Approaches" by Dan Simon, 2006 (ISBN: 978-0-471-70858-2).


The Sinc Function and the Partition of Unity Property

For any finite-energy signal $x(t)$, the function
$$y(t) = \sum_{n=-\infty}^\infty x(t-nT), -\infty < t < \infty$$
is a periodic function of $t$ with period $T$, that is, $y(t-T) = y(t), -\infty < t < \infty$. Note that $T$ is not necessarily the least period of $y(t)$.  Now, a periodic
function can be expressed as a Fourier series
$$y(t) = \sum_{m=-\infty}^\infty c_m \exp(j2\pi mt/T)$$ where
$$\begin{align}
c_m &= \frac{1}{T}\int_0^T y(t)\exp(-j2\pi mt/T)\,\mathrm dt\\
&= \frac{1}{T}\int_0^T \sum_{n=-\infty}^\infty x(t-nT)\exp(-j2\pi mt/T)\,\mathrm dt\\
&= \frac{1}{T}\sum_{n=-\infty}^\infty \int_0^T  x(t-nT)\exp(-j2\pi mt/T)\,\mathrm dt\\
&= \frac{1}{T}\sum_{n=-\infty}^\infty \int_{-nT}^{-(n-1)T}  x(s)\exp(-j2\pi ms/T)\,\mathrm ds\\
&=  \frac{1}{T}\int_{-\infty}^{\infty}  x(s)\exp(-j2\pi ms/T)\,\mathrm ds\\
&= \frac{1}{T}X\left(\frac{m}{T}\right)
\end{align}$$
where $X(f)$ is the Fourier transform of $x(t)$.
Now, for the case $T = 1$,
if $y(t) = 1$ for all $t$ as the OP desires, then it must be that $c_0 = 1$
and $c_m = 0$ for all $m\neq 0$, and therefore the Fourier transform of $x(t)$
must have the property that
$$X(m) = \begin{cases} 1, & m = 0,\\ 0, &m \neq 0,\end{cases}\tag{1}$$
that is, $X(f)$ is what is sometimes referred to as a Nyquist pulse (in the
frequency domain). We are familiar with Nyquist pulses in the time domain,
with $\operatorname{sinc}(t) = \frac{\sin(\pi t)}{\pi t}$ being the most
familiar one,  but there are many others, including as a special case
$[\operatorname{sinc}(t)]^k$. Applying duality,
$X(f) = [\operatorname{sinc}(f)]^k$ is a Nyquist
pulse in the frequency domain, and its inverse Fourier transform $x(t)$ is
a signal with the property that
$$\sum_{n=-\infty}^\infty x(t-nT) = 1, -\infty < t < \infty.$$
Note that $x(t)$ is the $k$-fold convolution of $\operatorname{rect}(t)$
with itself. For example, with $k = 2$,
$$x(t) = \begin{cases}1 - |t|, &-1 < t < 1,\\
0, & \text{otherwise.}\end{cases}$$
What about the other way around?  Well, $\operatorname{sinc}(t)$ has
Fourier transform $\operatorname{rect}(f)$ which has the Nyquist pulse
property $(1)$ in the frequency domain as does $[\operatorname{sinc}(t)]^2$
whose Fourier transform
$$X(f) = \begin{cases}1 - |f|, &-1 < f < 1,\\
0, & \text{otherwise,}\end{cases}$$
also satisfies $(1)$.  However, for $k > 2$,
the Fourier transform of $[\operatorname{sinc}(t)]^k$ does not
have the Nyquist pulse property $(1)$ because the three-fold (or more
generally $k$-fold convolution of $\operatorname{rect}(f)$ with itself
does not satisfy $(1)$: the $k$-fold convolution "bulges" outside
the interval $(-1,1)$.  See for, example, the
Irwin-Hall distribution
for details of the exact values of the convolution results.


Yes, the confusion comes from sampling theory. Thus, when you sample a signal, you get ... samples. If you want to be clear you can call the other one a statistical sample.


How does chessboard configuration affect camera calibration?

Generally speaking, when you are estimating a model from data, the more (independent) data points, the better the estimation.

In this particular case, the goal is to get as many 3d -> 2d point correspondences as possible, filling as much as possible of the volume of space of interest that is observed by the camera. The size of the squares directly affects how many data points you get per image (larger squares implies less visible corners/points).

It is also important to ensure that the camera parameters you want to calibrate are observable given the data. For example, in order to calibrate the focal length (or the field of view, equivalently) you need to observe a significant perspective foreshortening of the (planar) checkerboard in a few images. In order to accurately calibrate the nonlinear (e.g. barrel) lens distortion, you'll need to observe plenty of aligned data points near the edge of the image, where the distortion is larger.

There is a longer response on StackOverflow that you may want to look up.


How to model transfer functions using Python (without using higher level packages) and plot the magnitude and phase response?

If the response is real you need not convert it to complex in order to compute the magnitude, but the process is correct for complex and real number since the abs function will use complex conjugate product prior to taking the square root.  If it was real, there would be a trivial phase response since the phase would always be 0 degrees. This would not be a causal system so not a realistic example to have a real frequency response. For complex cases the approach used is correct for magnitude and phase, just not the most efficient.
To confirm results (or to compute directly), use the freqs and freqz commands that are part of scipy.signal for continuous-domain (s) and discrete-domain (z) transfer functions.
For example the transfer function
$$H(s) = \frac{s+1}{s+4s+4}$$ could be plotted using:
w,h = scipy.signal.freqs([1, 1], [1, 4, 4])

Returning the frequency (w) and complex response (h)
Read the help for both of those functions (also available under same name in MATLAB and Octave) for further info.
For MIMO systems the control package has this all worked out as part of the bode function.


Is my concept of receive diversity right? EGC,MRC and SC

Sounds pretty solid to me.

It's not necessarily MRC>EGC>SC, you'd need to write $\geq$ technically, since for the special case where all $|h_k|$ are equal you have MRC=EGC and for the special case of a single branch (antenna), they are all equal.


Why does cel-shading sometimes look similar to rotoscoping?

I checked out the Wikipedia page for rotoscoping and it says one way it was done is that a live action image was shown behind frosted glass and an animator would draw it.  The frosted glass detail obfuscation coupled with the fact that it was a cartoonist means to me that you are going to get solid colored areas as shading, like you'd see in a cartoon.  Cell shading is basically trying to make the same effect but using 3d rendered graphics as a source instead of an actual photographed image.  Seems like they are the same effect, just with different names!


The best way to get rid of 50Hz AC noise is with an Adaptive filter (usually the LMS algorithm will work fine), but you need to be able to take a reference signal of only the AC noise source, normally this is usually possible because if your system is being contaminated with AC it means you have some electrical equipment nearby or an outlet where you can take a sample.

If you do not want adaptive filter, then you can try a notch filter at 50Hz. I believe what your friend was trying to do is put a lot of zeros in the Region of 45-65Hz. It depends on the frequency response of the filter but normally the Harmonics will not be attenuated to the same degre, and normally AC is not very clean especially if you have transformers in the middle and such.

Adaptive filter is the way :)


Estimating the Variance of Noise Using Median Function Applies on the Gradient of the Image

This method is called median absolute deviation (MAD). Unlike the version with mean, it is a robust estimator of the variance. For a univariate set $X$, MAD is defined as:

$MAD = median_{i} (| X_i-median_j(X_j) |)$.

This is typically thought as the median of the absolute values of the deviations from the median of the data. The variance is then robustly estimated as:

$\sigma= \kappa MAD(X)$, where $\kappa=1.4826$ is a scale factor. Note that your $0.6745=\frac{1}{1.4826}$ corresponds to such factor.

In your example, MAD of the gradients are utilized to estimate variance of the noise. It is a pretty standard procedure. For more on the topic, please refer here.

In MATLAB, you could simply compute this via the function mad.


Reflecting negative frequency to positive frequency

You could use another approach: allow negative frequencies, do a complex IFFT, discard the imaginary parts of the time domain samples, and multiply the result by 2.
Let's try it in Octave (MATLAB clone) with a motif $[1 + 2i, 3 + 4i, 5 + 6i, 7 + 8i]$ shifted so that its leftmost bin ($1 + 2i$) lands on a negative frequency. FFT length is 8. (I rewrote the results for readability.)

> g = [3 + 4i, 5 + 6i, 7 + 8i, 0, 0, 0, 0, 1 + 2i]
g =

   3 + 4i   5 + 6i   7 + 8i   0 + 0i   0 + 0i   0 + 0i   0 + 0i   1 + 2i

> y = real(ifft(g))*2
y =

   4.00000  -0.89645  -2.00000   0.98223   1.00000  -1.60355   0.00000   4.51777

> fft(y)
ans =

   6 + 0i   6 + 4i   7 + 8i   0 + 0i   0 + 0i   0 + 0i   7 - 8i   6 - 4i

The same time-domain signal $y$ can be generated by synthesizing the first 5 bins of $\text{fft}(y)$ and feeding them to real IFFT (which I assume you use). The 3 remaining negative-frequency bins the IFFT kind of generates internally. You can follow these steps for the synthesis:

1) For each bin of the shifted motif that lands on a positive-frequency bin, copy the value verbatim. (You are already doing this)

2) For a bin of the shifted motif that lands on the 0 Hz bin, double its real part and write that to the 0 Hz bin, discarding the imaginary part, like $3 + 4i \rightarrow 6 + 0i$ in the example. (You are already doing this)

3) For each bin of the shifted motif that lands on a negative-frequency bin, add its complex conjugate (same real part, sign of imaginary part flipped) to the positive-frequency bin that is at the same distance from the 0 Hz bin. This will be one of the bins already written to in step 1. In the second bin of the example $5+6i$ and $\text{conj}(1 + 2i) = 1 - 2i$ were summed to get $6 + 4i$. Don't expect the imaginary parts to cancel, because even if the motif might be symmetrical, those two motif bins are not located symmetrically on the opposite sides of the motif, unless you are synthesizing a 0 Hz sinusoid. It is perfectly fine if the imaginary part remains non-zero; frequency domain bins can have complex values and still yield real-valued time-domain data.

Oh, and if the shifted motif reaches the Nyquist bin (corresponding to half the sampling frequency) or beyond, do exactly the same things there, mirroring around the Nyquist bin.


non-consistent normalization problem in fft

I believe you need to scale by the sample period $\Delta T$ properly in your differentiator (divide by $\Delta T$) if your time units are in seconds not samples.

Your differentiator can be arrived at by using the backward Euler approximation of a continuous time differentiation.

The mapping from s to z using the backward Euler

$$s \leftrightarrow \frac{1-z^{-1}}{\Delta T}$$

$1-z^{-1}$ is simply subtracting the past sample from the current sample. To maintain your time axis in seconds, divide by $\Delta T$. If your time axis is in samples (as we often do), the $\Delta T=1$ and further scaling is not needed.

Similarly an integrator which is $1/s$ would be multiplied by $\Delta T$. For all other mappings you see the same $\Delta T$ appear, just as it exists as $dt$ in the continuous time equations.


MATLAB: Adding Noise with Regard to Signal to Noise Ratio (SNR) or EbNo?

These stuffs $E_s/N_0, E_b/N_0 \textrm{ and SNR}$ are convertible.

\begin{align}
E_s/N_0 &= E_b/N_0 + 10\log_{10}(k) \\
E_s/N_0 &= 10\log_{10}(T_{sym}/T_{samp}) + \mathrm{SNR}
\end{align}

where $k$ is the number of information bits per symbol, $T_{sym}$ is the signal's symbol period and $T_{samp}$ is the signal's sampling period. More details can be found at AWGN channel MATLAB.

And by description of awgn() function:

  y = awgn(x,snr) adds white Gaussian noise to the vector signal x. The
  scalar snr specifies the signal-to-noise ratio per sample, in dB. If x
  is complex, awgn adds complex noise.

In your case, no oversampling thus awgn() uses its snr parameter to generate $E_s/N_0 =$ snr.

You are free to plot your BER curves over $E_s/N_0$ or $E_b/N_0$ because they are equivalent. You just need to convert them back and forth via the modulation order $k$.

However $E_b/N_0$ is usually prefered because (EbN0 wikipedia)

  $E_b/N_0$ directly indicates the power efficiency of the system without
  regard to modulation type, error correction coding or signal bandwidth
  (including any use of spread spectrum). This also avoids any confusion
  as to which of several definitions of "bandwidth" to apply to the
  signal.


The fft() function calculates the one-dimensional FFT of its input argument. If the input argument is a vector, then the operation is pretty simple to understand; the output is just the result of efficiently calculating a discrete Fourier transform on the input. If the input is a matrix, then, as many MATLAB functions do, each column is treated as a separate one-dimensional signal and are transformed separately. Therefore, the $n$-th column of the output matrix corresponds to the DFT of the $n$-th column of the input matrix.

The fft2() function calculates the two-dimensional FFT of a matrix. One way of envisioning this operation is to take standard one-dimensional FFTs down each column of the input matrix $A$ to yield a temporary result $B$. Then, take one-dimensional FFTs across each row of $B$ to yield the result $C$. The matrix $C$ would then be the result of the two-dimensional FFT of A. 2-D transforms (although not typically the FFT) are often used in image and video compression algorithms.


The classic complexity measures are much less relevant today than they were in the days of fixed point ALU's taking 40 cycles to do a floating point multiply.

The answer to your question is maybe.  One can decompose an FFT in terms of smaller FFTs,

based on

  Rabiner, Lawrence R., and Bernard Gold. "Theory and application of
  digital signal processing." Englewood Cliffs, NJ, Prentice-Hall, Inc.,
  1975. 777 p. (1975).

there is a section that shows how to do a 1 dimensional dft in two dimensions.

see my answer at

Combine FFT's of shorter length than sample data to get spectrum of all data

and coppying the matlab code from there:

for a 96 point FFT, I formed a 32 X 3 matrix

clear all
M=3;
N=32;
x=linspace(0,10,M*N);
X=reshape(x,N,M).'; % read in as rows
Twiddle=zeros(size(X));
for i=1:M
for k=1:N
Twiddle(i,k)=exp(-1j*2*pi*(i-1)(k-1)/(NM));
end
end
X=fft(X) % fft on each column
X=X.*Twiddle;% element by element product
X=fft(X.').' ; %fft on each row
y=reshape(X,N*M,1); % read out as columns
figure(1)
plot(abs(y),'linewidth',2)
title('Composite DFT')
figure(2)
plot(abs(fft(x)),'linewidth',2)
title('Direct DFT')
figure(3)
plot(x,'linewidth',2)
title('time series')

You need to exploit the zero padding when you form the 2D Matrix. When you do the column FFT operation, simplify by noting the zero elements.

As a simple example, let's say we have 3 data points we want to zero pad for a 16 point DFT. A direct DFT is complexity $16^2$. We need to choose M and N, so that M*N=16. I want all the non zero data in the first row, so M=N=4 , so read in as rows is:
$$
X=\left| \begin{array}{cccc} x[0] & x[1] & x[2] & 0 \\
 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 \\
\end{array}\right|
$$
the next step is to do DFTs on the column, normally $ 4 \times 4^2$ complexity. The last column are zeros, and the DFT of a column is zero, so
our complexity is now $ 3\times 4^2$. Im reusing the same memory,
$$
\text{after column DFTs} \; \;X=\left| \begin{array}{cccc}
 x[0] & x[1] & x[2] & 0 \\
 x[0] & x[1] & x[2] & 0 \\
 x[0] & x[1] & x[2] & 0 \\
 x[0] & x[1] & x[2] & 0 \\
\end{array}\right|
$$
There was only one term in each column so I didn't really need to do any multiplies in this step and you don't either.
The Twiddle matrix is usually calculated offline in any FFT , and  forming it is not usually counted in the Complexity calculation.
In most DSP treatments, the top of the array is zero based, but in Matlab is one based.
so in MATLAB
$$ W_{16}^{ik}=
\exp(-\jmath 2\pi\frac{(i-1)(k-1)}{16})
$$
but in most DSP books
$$ W_{16}^{ik}=
\exp(-\jmath 2\pi\frac{(i)(k)}{16})
$$

In this case the product of the Twiddle matrix with $X$ is:
$$
\text{Twiddled X} X=\left| \begin{array}{cccc}
x[0] & x[1]1 & x[2]1 & 0 \\
x[0] & x[1]W_{16}^1 & x[2]W_{16}^2 & 0\\
x[0] & x[1]W_{16}^2 & x[2] W_{16}^4 & 0\\
x[0]& x[1]W_{16}^3 & x[2]W_{16}^6 & 0\\
\end{array}\right| \;
$$
given
$$
\text{Twiddle}\; X=\left| \begin{array}{cccc}
1 & 1 & 1 & 1 \\
1 & W_{16}^1 & W_{16}^2 & W_{16}^3 \\
1 & W_{16}^2 & W_{16}^4 & W_{16}^6\\
1 & W_{16}^3 & W_{16}^6 & W_{16}^{9} \\
\end{array}\right| \;
$$
at this step 9 Complex products

The rest should be obvious (and is left as an exercise)
Use the Matlab code as a template but do the short cuts when you can.


Circular buffer wrap around with chorus effect

If I understood correctly what you are asking, then this is probably more of a programming question... But the solution would be to use modular arithmetic to compute the position of the delay. With a length equal to a power of 2, this could be accomplished with an AND instruction.

const int Mask = (1 << 16) - 1;
Circular_Array[(Current_Position_in_Buffer - Delay) & Mask];

In this example, the length of your buffer is $2^{16}$.

More generally you can use the modulo operator, but it will be slower.


Scilab's filter, for short coefficient vectors, function implements a linear convolution in C code; that alone, since there's no python to actually be evaluated here, just multiplication and addition, is much much faster than writing something in a scripting language that can't 100% be just-in-time compiled.

For longer vectors, scilab implements fast convolution; ie. it exploits the fact that (circular) convolution in time domain corresponds to point-wise multiplication in (discrete) frequency domain, and uses zero-padding and saving of overlaps to emulate the linear convolution (which filtering represents) with that.

So, either way, use your libraries when doing signal processing! Aside from the convolution, there's other things that are generally faster if done via clever usage of library functionality: For example, whenever you have a loop that looks like

sum = 0
for a, b in zip(vectorA, vectorB):
    sum += a*b

you'd be far, far better of doing a dot product of the two vectors.

You have to consider this: Your CPU is very fast at doing basic math operations – often, it can do for example 8 multiply-and-accumulates (MAC) operations in a single step. Compared to that, parsing the structure of the (precompiled, even) for loop, building temporary python objects to hold the individual values for a and b, and overwriting the sum object, thus removing the old object and replacing it with a new one, leading to garbage collection and so on, is way way way way more work than just doing the maths. I like to put it like this:

  Imagine you're tasked with multiplying a lot of numbers between 0 and 10, but the numbers you need to multiply are written in text form in a book.
  Reading that book will take much, much longer than the multiplications

That's how it is to use dynamic languages to do basic math operations.


You have a very narrow stop band which means that all the poles are crammed in a very small area of the complex plane, close to the unit circle. This can result in severe numerical problems, even for relatively small filter orders, even with floating point arithmetic.

Another important point that you might not realize is that if you design a band pass or a band stop filter using the command cheby2 (and all other similar commands), then N is the order of the prototype low pass filter, and not the order of the resulting band pass or band stop filter. So if you choose N=14 then you've designed a band stop filter of order $28$. Such a filter is very hard, if not impossible to realize because of numerical problems.

Also note that calling cheby2 with [z,p,g] (zeros, poles, and gain) as output argument usually results in higher accuracy than using the polynomial coefficients.

I would actually suggest trying a notch filter instead of such a narrow band band stop filter.


AR Modeling: Why residual is white noise?

it's a hand-wavy argument, but the idea is to derive the coefficients $\alpha_i$ so that the norm of the error $\Big\||\epsilon[n]|^2\Big\|$ is minimized.  we assume that $x[n]$ is a stationary "random" process, we know all of the previous samples $x[n-i]$ and we want to make a good guess at $x[n]$ with a linear combination of the previous samples.  we can only make some guess of the next $x[n]$ if it is not a white random process.  if it is white, knowing the previous samples is not going to help us at all.  if it is not white, then the autocorrelation of $x[n]$ is not zero for non-zero lags and from that we can compute $\alpha_i$.

now if $\epsilon[n]$ was not white, we could do further LPC to make a good guess at $\epsilon[n]$ and make a better guess for $x[n]$, but if the coefficients $\epsilon[n]$ are optimal, then we can't make a better guess for $x[n]$.  so then the error of our guess for $x[n]$ and what we really get has to be white.


Detection of wobbling sound with rising notes

I would start with following set of parameters:

MFCC's (I that know you tried it, but stay with me) without static energy (1'st coefficient)
Some descriptors from MPEG-7, like: Spectral Flatness, Harmonicity, Fundamental frequency, Spectral Spread, etc. You can find out more about them here: click or in this great book: MPEG-7 Audio and Beyond

You should also include $\Delta$'s and $\Delta\Delta$'s of MFCC's and Spectral Flatness. They should really improve accuracy. You can also try to calculate these for other parameters. Having that, you can start with Gaussian Mixture Model or Support Vector Machine (I think that GMM will be suitable for that task). After training your model simply validate it on testing set (seems obvious but I've seen many projects where people tested their Machine Learning algorithms on training set with added noise or, what's even worse, same training data).

After this, you should perform Feature Optimisation. This is important part, as some of the features are irrelevant (i.e. Spectral Flatness for few first filter banks, while bird's chirp has higher frequencies). Nice publications with many algorithms for that purpose can be found here: click,   click.
Good luck!


Simulating Range Bearing Sensor with MATLAB with Gaussian Noise (Generating Gaussian Colored Random Vector)

Few notes:

First you should multiply the noise by the Standard Deviation (Root of the Variance for zero mean noise).
You can do that by multiplying the Lower Cholesky Decomposition of matrix by a column vector of Gaussian noise. Yet since you assume no correlation, you can do that by independent multiplication.

Something like:

vZ = vZ + [sigma_range; sigma_angle] .* randn([2, 1]);


Coherent Sampling vs Windowing for 50/60 applications

That's coherent sampling.

Assuming this is a AC power application: it's only coherent sampling if you derive the sampling clock directly from the AC frequency using a PLL (Phase Locked Loop) or equivalent.

But does this account for frequency variations?

Depends on how you derive the sampling clock. If you use a PLL it will track the AC frequency, if you use an independent clock, the frequency will indeed drift.

Should I use windowing? What other approaches can be implemented?

If you really just want to measure the RMS of the signal, you don't need coherent sampling. Actually incoherent is better. You just need to make sure you have enough samples.


Create a color image in Matlab with sinusoidal patterns

Simply do it for each color channel, I mean for Red,Green and blue, for each form a separate matrix and then concatenate those, to form a 3D color image.

fs = 0.08;   //FOR RED CHANNEL
W = 256;
H = 128;

img1 = zeros(128, 256);
for m = 0 : H-1
    for n = 0 : W-1
        img1(m + 1, n + 1) = sin(2*pi*fs*m);
    end
end

fs = 0.16;    //FOR Green CHANNEL

img2 = zeros(128, 256);
for m = 0 : H-1
    for n = 0 : W-1
        img2(m + 1, n + 1) = sin(2*pi*fs*m);
    end
end

fs = 0.01;    //FOR Blue CHANNEL

img3 = zeros(128, 256);
for m = 0 : H-1
    for n = 0 : W-1
        img3(m + 1, n + 1) = sin(2*pi*fs*m);
    end
end

Color=zeros(128,256,3);
Color(:,:,1)=img1;
Color(:,:,2)=img2;
Color(:,:,3)=img3;
imshow(Color,[]);


Signal 1 was an average of a left and right signal and signal 2 was a binocular signal.

I followed the suggestion of MBaz, and computed the frequency response of the filter for each segment (N = 3) for each subject (N=12).  Then I averaged the frequency responses of the filters. Here is the result:

MBaz and Matt: Thank you so much for solving my problem.
Matt, it looks close to your filter.
Lee


A signal can be decomposed into the sum of an even (symmetric) signal plus an odd (anti-symmetric) signal.  One measure of symmetry of the signal might be the ratio between of the energy in the even decomposition to the energy of the total signal.


After much searching and asking questions, I have written my own realization for the vibrato effect: https://github.com/Bershov/Vibrato-effect. The code is not hard for understanding and works fine. In fact, the vibrato is a kind of frequency modulation. Its key concepts: variable delay, low-frequency oscillation and intelpolation.


The ratio of sample playback rates should be equal to the ratio of pitches you want to obtain.

For example, if your C4 note was sampled at 48kHz, you'll need to play it back at $48000 \times \frac{277.18}{261.63} = 50.85kHz$ to make it sound like a C#4.


What factors affect the propagation delay at wireless communication?

Even if transmitter and receiver are fixed, you have unknown reflecting objects, resulting in unknown multi-path propagation. The resulting propagation delay can be modeled as a random variable. If there are moving objects in the environment (such as people moving in a room), then the propagation is not only unknown but also time-varying. The delay spread of a channel describes the difference between the maximum and the minimum propagation delay in the given environment.


Meaning of channel capacity in the AWGN channel

There are two different kinds of channel models that are being confused here.

In the binary symmetric channel, the inputs and the outputs
are constrained to be $0$ or $1$ and the key parameter is the
transition probability: the probability that an input $0$ is changed
to an output $1$ or vice versa.  The channel capacity $C$
is a number between $0$ and $1$ and the meaning is that by using
a suitable code of rate $R < C$ on the channel, we can transmit
at a rate of $R$ information bits per channel use with arbitrarily
high reliability (that is, with arbitrarily small bit error
probability.  A (block) code for such a channel would be a set
of $M$ binary $n$-tuples that will convey $\log_2 M$ bits in $n$ channel
uses at a rate $R = \frac{\log_2 M}{n}$. Usually $M$ is chosen as
a power of $2$, say $2^k$ for some $k < n$.  Note also that the
binary symmetric channel capacity cannot exceed 1 bit per channel use.

The discrete-time Gaussian channel is a different beast. The inputs
can be arbitrary real numbers (it is not necessary that they
be restricted to $0$ or $1$) and if the $i$-th input is $x_i$,
then the $i$-th output is $x_i+n_i$ where the $n_i$ are independent
zero-mean Gaussian random variables with variance $\sigma^2$.
$x_i^2$ is called the power in the $i$-th input and
$\sigma^2$ is called the noise power.  A (block) code for such
a channel is a collection of $M$ $n$-vectors $\mathbf x^{(i)}$
(called a codebook) that can be used
to transmit $\log_2 M$ bits in $n$ channel uses.  Various realistic
constraints are placed on the code.  One
such constraint is that the peak power is bounded by $P$:
the magnitudes of all the $M\times n$ entries $x_j^{(i)}$
in the codebook are no larger than $\sqrt{P}$.
Another is that
that the average power is bounded by $P$:
$$\frac{1}{Mn}\sum_{i} ||\mathbf x^{(i)}||^2
= \frac{1}{Mn} \sum_i \sum_j \left|x_j^{(i)}\right|^2\leq P.$$
The capacity of the channel depends on the signal to noise ratio
$P/\sigma^2$ and it can exceed one bit per channel use.  Think of it
this way: we can use multilevel pulse amplitude modulation to transmit
more than one bit per channel symbol on this channel whereas with
the binary symmetric channel we are restricted to binary modulation
and so cannot exceed one bit per channel use.

For lots of details on this kind of stuff at a level lower than
a typical graduate text, see the following links (lecture notes
of mine).

Lecture 6,
Lecture 7,
Lecture 8,
Lecture 9,
Lecture 10,
Lecture 11,
Lecture 12


Fourier transform of even/odd parts of a complex signal

Defining
$$ x_e(t) = \frac{x(t) + x(-t)}{2} $$
$$ x_o(t) = \frac{x(t) - x(-t)}{2} $$
We can take the Fourier transform of each:
$$
FT\{x_e(t)\} = \int_{-\infty}^{+\infty} x_e(t) (\cos(2\pi\omega t) - j \sin(2\pi\omega t)) dt
$$
$$
FT\{x_o(t)\} = \int_{-\infty}^{+\infty} x_o(t) (\cos(2\pi\omega t) - j \sin(2\pi\omega t)) dt
$$

So if $x_e(t) : \mathbb{R} \mapsto \mathbb{R} $ then the $\sin$ term disappears and we are left with just a real-valued function.

So if $x_o(t) : \mathbb{R} \mapsto \mathbb{R} $ then the $\cos$ term disappears and we are left with a purely imaginary-valued function.

Now, if instead $x_e(t), x_o(t) : \mathbb{R} \mapsto \mathbb{C} $ then the term cancellations still happen, but we are still left with a complex value for each integral because both functions are now complex-valued.

As a result, the right-hand sides of your quoted equations do not hold in general.


Note that for stable IIR filters, the impulse response does approach zero as $n$ goes to infinity. It just never becomes exactly zero. However, the sum of the absolute values is finite. Just as an example, take the exponential impulse response

$$h[n]=a^nu[n],\qquad |a|<1\tag{1}$$

where $u[n]$ is the unit step function. The sum

$$\sum_{n=-\infty}^{\infty}|h[n]|=\sum_{n=0}^{\infty}|a|^n=\frac{1}{1-|a|}\tag{2}$$

is obviously finite (for $|a|<1$).

I agree with Hilmar's answer, but I would like to add that everything he said about poles is only valid for causal filters. In the general case it is not necessary for a stable filter to have all its poles inside the unit circle. For stability it is necessary and sufficient that the unit circle is inside the region of convergence (ROC) of the filter's $\mathcal{Z}$-transform (transfer function). So a stable filter can have poles inside as well as outside the unit circle, but in that case it can't be causal. So it's only the additional requirement of causality that forces the poles of a stable IIR filter inside the unit circle.

A similar thing holds for FIR filters. Causal FIR filters indeed have all their poles at the origin, but non-causal FIR filters also have poles at infinity. A simple example is this non-causal moving average filter

$$H(z)=\frac15\left(z^{-2}+z^{-1}+1+z+z^2\right)\tag{3}$$

which has poles at $z=0$ as well as poles at $z=\infty$, but being an FIR filter, it is obviously stable.


Add harmonics to a FFT decomposed signal

For each added harmonic, for each bin number below N/2, you could divide by the harmonic ratio, then use that resulting number to interpolate a point below in the FFT, unless the division results in an integer, in which case you could just use the value from a lower FFT bin.  Use a high-quality interpolator, such as a windowed Sinc or other polyphase upsampling filter kernel.  For a strictly real result from an IFFT, make sure the values above N/2 in your harmonic frequency synthesis remain conjugate symmetric with those below N/2.


Time signal FFT different for whole signal than half signal on high frequencies

The primary issue is due to using a rectangular window with the floor being the result of spectral leakage.

This is demonstrated in the plot below showing the frequency response (kernel) for the two rectangular windows used; one 500 samples long as red and the other 1000 samples long as blue. The corresponding rectangular window is effectively multiplied by your signal in time, thus convolves with your signal in frequency, resulting in spectral leakage as the spectrum given in the original question.  Thus you are seeing the smearing (convolution) of the actual spectrum with this kernal in each case. Using a window with a wider dynamic range (lower sidelobes) such as the Blackman  would result in a better view of the actual spectrum to the extent that spectral components now are getting buried under the spectral leakage.

For further details on windowing and window choices, please see this post: Why would one use a Hann or Bartlett window?

You will also see a 3 dB difference between two windows if the signal energy is distributed wider than the resolution bandwidth in each case, since the resolution bandwidth for the rectangular window is $1/T$ where T is the length of your waveform. (So when you double the length from 500 to 1000, the resolution bandwidth is half resulting in half the power in each bin (unless a spectral tone occupies one bin in which case you would not see a change).


This would be an allpass filter. Except for the trivial case of unity and integer-sample delays, these can't be done as FIR filters and in general an IIR filter is required. However, they are easy to make. The zeroes of an allpass are simply the inverse of the poles (and vice versa). If you have the poles in polynomial form, you can simply flip them to get the zero polynomials. For example a second order allpass looks like this
$$H(z)=\frac{a_{2}\cdot z^{0} + a_{1}\cdot z^{-1} + a_{0}\cdot z^{-2}}{a_{0}\cdot z^{0} + a_{1}\cdot z^{-1} + a_{2}\cdot z^{-2}}$$
Strict allpass filter have $\left \|H(e^{j\cdot \omega })\right \|=1 $ for all frequencies. You can certainly design approximation using FIR filters if you only need this property for a limited frequency range and if the magnitude doesn't have to be exactly unity.


CIC Filter implementation for Normalized Correlation

So, assuming you already know that the two streams you're multiplying point-wise are aligned, and you want the output of your system to be the correlation coefficient of zero shift for the length of the CIC-built averager, this is the way to go.

So, the output of your CIC will be the sum of all point-wise products of the last N samples flying through.

Now, if either input stream becomes overly "strong", that product might become very large, regardless of whether the correlation is actually strong or not.

To avoid that, you'd usually just divide by the square root of the product of energies of both sequences. Assuming the signals are stochastic stationary signals, that would be the product of the square roots of N times the variance, which is simply N times the product of the standard deviations. If you have knowledge of these, a simple multiplication of an input or the output by 1/(N·std1·std2) would do.

Assuming you have no knowledge, you can measure: Just use the same CIC structure on the two squared inputs (instead of their product) to get the instantaneous energy, and divide the result of your CIC correlator by the square root of that.


Sum of sines isn't a concrete waveform, there is no general form of its ambiguity function. But for any given sum of sines you can calculate this function by evaluating cross correlation for various time and Doppler deviations. Though one rule exists, if the number of sinewaves in composition grows, the ambiguity function will become narrower in both frequency and time dimensions, since such a signal will be more and more noise like.

Since autocorrelation isn't linear function of some signal $s$ (its quadratic as you know cause it is sum of multiplications of $s$ and its shifted copy), so you can't express autocorrelation and more generally ambiguity function of sines composition as a sum of their autocorrelations (or amb. func-s.). Actually it is the sum of entering sines' autocorrelation functions and all possible entering sines' cross correlation functions. It could be something huge. So the most general way is to compute of such a function for given sines sum directly with knowing formula of autocorrelation.


position of spike in the spectrum of a sin wave

A sinusoid is the sum of two complex exponentials ("negative" and "positive" frequency) so that their imaginary parts cancel out each other:

$$cos(f \cdot t) = \frac{1}{2} \left( exp(f \cdot t \cdot i) + exp(-f \cdot t \cdot i) \right)$$

The Fourier Transform decomposes your input into a weighted sum of complex exponentials. Therefore, you get to see two pulses in the spectrum. One at -f and another at +f, basically. Typically, the "negative half" of the spectrum is located after the positive half because in discrete signals there is no difference between -f Hz and fs-f Hz where fs is the sampling frequency. So, your FFT result's ordering is 0 Hz to almost fs Hz where you can think of the second half as the negative portion of the spectrum.

Since real-valued signals have an imaginary part of zero, you'll notice that the Fourier Transform of real-valued signals always exhibits a symmetry between negative and positive frequencies with a certain kind of phase relationship (same amplitude and negated phase) so that the complex exponentials' imaginary parts cancel between negative and positive frequency.

Complex-valued signals are more general in that their spectra don't have to be symmetric around 0 Hz. The imaginary part helps you distinguish between positive and negative frequencies:

$$exp(f \cdot t \cdot i) = cos(f \cdot t) + i \cdot sin(f \cdot t)$$

As you can see above, the sign of f does not change the real part. But it does affect the sign of the imaginary part.


Advantages and disadvantages of Hilbert-Huang Transform versus Fast Fourier Transform?

Fourier algorithm is actually a global transform that can not reflect the damping and local specialty. Unfortunately, these features are fundamental and important for characterize a nonstationary signal. on the other hand in the Fourier algorithm the estimation of frequencies are sensitive to noise. so HHT is better in your case if you are going to deal with noisy signals.


Tetrahedral microphone array beamforming

I don't know what the fuss is about.  This is a simple vector problem.

"Well, if it is that easy please elaborate. I have 4 microphones on the following locations (0.042,45,35),(0.042,-45,-35),(0.042,135,-35),(0.042,-1355,35). These are coordinates for the microphones with m1=(ri,θi,φi), where −π≤θ≤π and −π2≤φ≤π2. create a beam towards a sourse at (θs,φs)=(22.5,67.5). – havakok"

You need to define five vectors in Cartesian coordinates, one for each mic and one for the beam.  The beam vector should be normalized to a length of one.

If the sound source is sufficiently far away, the assumption that the wavefront is coming in as a planar wave is a fairly accurate approximation.  Therefore, you simply have to find the projection of each mic location onto the line defined by the beam vector.  The projection for each mic is found by dotting the mic vector with the normalized beam vector.  The value you get is the distance along the line (where one is the beam vector length).  Use the speed of sound to convert this distance to a time shift.  Use your sampling rate to convert the time shift to a sample shift.  Shift and add.  The last mic to receive the signal is not delayed, and the previous mics are delayed relative to their projection's distance to the last mic's projection.

It should be obvious that you only get one beam direction at a time.  Don't be oversold on the "beam" notion either.

Ced

Two Answers:

1) Rounding to the nearest sample should be negligible.

Some ball park figures to consider:  Suppose you have a 1000Hz tone.  Sound travels at roughly 1130 feet per second.  Thus one cycle of the tone is roughly 1.130 feet long, or about 13 and 1/2 inches.  At CD quality sampling, about 44 samples.  So each sample is about 3/10 of an inch (this is frequency independent).  Compared to the length of the cycle, this isn't a lot.  However, if you are concerned about frequencies above 10k, this may matter.  In that case, a simple linear interpolation in your delay should do the trick.

2) The beam effect is formed by "coordinated interference effects", i.e. signals coming from other directions will be added together in an out of phase manner (only the signals in the desired direction are in phase), thus you are relying on interference to destructively interfere.  This is frequency dependent.  Also, signals coming in nearly at the beam angle won't be affected much, so perhaps "cone forming" is more accurate.

You didn't mention whether these were omnidirectional mics or not.  No bearing on the math, but may influence your results.


Terminologies - sparse channel, sparse input, compressed sensing

The term sparse, as you mention, refers to the fact that some "signal", usually represented by a vector $x$ contains mostly zero or negligible values and only a few non-zero or significant values. In fact, the situation with neglible vs significant values which is not strictly sparse is often referred to as "compressible", see e.g. Candès & Wakin 2008 (you can easily search for an openly accessible PDF of that paper).

Without knowing the literature you are referring to I can only make some guesses, but in general this is what I would assume:

The channel being sparse usually means that its (discrete) impulse response is sparse. In the case of speech recording, my guess is that they assume the impulse response from the speech source to the microphone can be modelled as a few distinct reflections of the sound off the surrounding walls.
A signal in general is sparse whenever you can find a sparse representation of it in some dictionary. That is, if we have some signal $x$, it is not necessarily $x$ itself that has to consist of only a few non-zero coefficients and a large number of zeros. Whenever we can find some "dictionary" $D$ such that $$x = D z$$ where $z$ is a sparse vector, we can say that $x$ has a sparse representation and compressed sensing can in principle work.
The term compressed sensing, very generally, can be used whenever you are trying to solve an under-determined system of linear equations $$y = M x = M D z = A z,$$ where $A$ is $M \times N$, where $M<N$ and you exploit the fact that you can actually solve it because you know that the number $K$ of non-zero coefficients in $z$ is $K<M$ (simplified - $M$ has to be $O(K \log N)$. It can be applicable to input and channel alike, depending on which one is assumed sparse. The signal model used for reconstruction of the sparse quantity just has to be posed accordingly, whether it is the speech signal or the channel through which it is recorded that is assumed sparse.

It seems you are being confused by the fact you have most often seen compressed sensing applied to the recorded signal so that $x$ above is the speech recorded by the microphone and the impulse response would then be represented by the measurement matrix $M$ above or, using the variable $H$ - for impulse response matrix - instead: $$y = H x$$ In the case you mention it sounds as if the channel may be what they are interested in sensing instead. I that case one could imagine recording a known signal $x$ in order to estimate the channel (impulse responce), in which case we can swap the roles of channel and signal and instead use the following signal model: $$y = X h$$ to be able to estimate the (sparse) channel $h$ from measurements $y$ and our known speech signal $x$ used to compose the Hankel matrix $X$.


describing a block diagram // How do I describe multiplying a signal through multiple branches formally?

What I would do is:

Add $m$ multiplier blocks, using a circle with a $\times$ inside
Each multiplier would have inputs $X$ and $Y_k$
Each multiplier would have output $f_k$
Be consistent in the use of bold typeface, uppercase and lowercase

Then I would say something like "The system has input $X$ and outputs $f_k=XY_k$, for $k = 0,\ldots,m$. In the diagram, $\times$ stands for signal-matrix multiplication, which is defined as follows...".


There is not really a difference between wrapped and unwrapped phase, it is mainly more visually pleasing when it is continues. The only way this could occur in an actual transfer function is when you would have two complex conjugate unstable poles and two complex conjugate stable zeroes with the same complex and real part, but the real part with opposite sign and the real part should be very small.

One way you could unwrap the phase is by splitting the transfer function up into multiple terms. Since multiplying two transfer functions, means adding their magnitude in dB and phase in degrees of radians. So you could first calculate the bode plot of the two complex conjugate poles (which should stay within 180°, unless it already wraps -180° to 180°) and then add the bode of the real pole to it. Or in this case you could look when the phase is bitter than 0°, in which case you subtract 360°.

In order to improve you model I would take a slightly better at the pole location, since their phase does not seem to drop instantly, so a small real part. There also might be a small delay ($e^{-Ts}$), since the phase of the data drops further than -270°.


How to get Instantaneous Magnitude for a Instantaneous Frequency From FFT?

Frist of all you need compute the Magnitude:

windowed = framed .* hann(length(framed))
Fourier = fft(windowed)
Mag = abs(Fourier)

And you need build a Window Kernel based in your Window Function!

For a Hann Window the kernel can be:

   Wk(k) = 0.5 * (sinc(k*(M/N)) / (1 - (k*(M/N))^2))

M=Window Size
N=FFT Size

Of course you can use the same size for both, N and M !

k = Bin Number from FFT / N
sinc = Normalized Sinc, visit http://en.wikipedia.org/wiki/Sinc_function

If you have already determined the Instantaneous Frequency, you need find the corresponding offsets for each Frequency and the Instantaneous Magnitude is found when you multiply the magnitude from FFT by your Hann Window Kernel index:

 Offset = abs(InstFreqs(i) / (Fs/N) - i)
 index = floor(abs(Offset) * N) + 1
 InstMag = Mag(i) * Wk(index)


Why does a longer observation time improve DFT resolution, but repeating a signal does not?

Why is this not equivalent to simply observing the signal for 1 period, and then paste it together N times?

It's only equivalent if certain conditions are met. Let's look at a single sine wave with frequency $f_0$ and period $T_0 = 1/f_0$. Sample rate is $f_s$, sampling period $T_s = 1/f_s$ and DFT length is $N$.
The input to the FFT would be
$$x[n] = sin(2\pi \cdot f_0\cdot T_s \cdot n), 0 \leq n < N$$
Now let's create two signals of length $2N$: one with double the length and one with pasting two buffers together.
$$x_1[n] = sin(2\pi \cdot f_0\cdot T_s \cdot n), 0 \leq n < 2N$$
and
$$x_1[n] = sin(2\pi \cdot f_0\cdot T_s \cdot mod(n,N)), 0 \leq n < 2N$$
The only difference is that modulo function in the second one, but it's important. It basically resets the counter to 0 when it reaches an integer multiple of N. In general this will create a discontinuity in the signal. It forces $x_2[N] = x_2[0] = 0$ even if $x[N]$ would be non zero. The first halves of $x_1$  and $x_2$ are the same but the second halves are not. $x_1$ is still a sine wave but $x_2$ isn't anymore.
The ONLY case where pasting is the same as making the signal longer is if the period of the signal is an integer multiple of the sampling period. In this case, there is no discontinuity and both cases are identical. And sure enough, in this case the DFT will resolve two different frequencies perfectly (as long as both meet the integer multiple condition).
A more philosophical interpretation: Making the signal longer provides more information. Pasting buffers just repeats information you already have, so you are not gaining anything new.


Why does part of my image turn black when applying functions to a part of an image in Matlab?

'edited image' is of type uint8 and the selected area which you have taken for processing finally after applying all operations is coming out as type double.But when you place it in original image it is considered as of type uint8 because the image is of type uint8. Hence before putting it in original image you should multiply x by 255 so that it get be properly scaled.


Cannot figure out the convolution operation in this problem

Hint:

For which values of $k$ is $x[k]=\delta[k-1]$ not equal to zero? How does that simplify the given sum?


I don't believe there is a right answer.  DCT or IDCT will achieve the same purpose: decorrelating to put most energy in fewer coefficients.  Whatever you do next with your MFCC (compression, feature extraction) will work with either approach.

Similarly, I think there are implementations that do an FFT instead of an IFFT for the first stage of the MFCC extraction.


How is the topo / scalp map from an EEG Signal with ICA created?

In general, the topographic map describes the distribution of electrical activity across the brain as measured in a number of well known "stations" (electrodes whose position is known) on the scalp.

This electrical activity is the combined result of millions of neurons firing below the skull. The strength of the electrical signal they create is measured across the REFerence electrode and one of the scalp electrodes. In a typical Electroencephalography (EEG) montage, you may have something like 16 to 20 channels (more in special cases).

The Inverse Problem in this case, that is, given the 20 (for example) signals of the EEG, calculate which neurons fired and how is ill-posed because it has a huge number of unknowns (millions of neurons) and only a handful of knowns (20 signals).

What can be done however, is to use Independent Component Analysis (ICA) to decompose the combined electrical activity of millions of neurons in the brain to a small number of independent (in a very specific mathematical sense) components.

In this case, the physical interpretation of this process is that the 20 channel activity is produced by 3,4,5,... independent "sources" within the skull. And, depending on the cognitive task that the subject was going through during data acquisition (whose performance is modulated by their mental state), these independent components might be able to provide some additional information. For example, maybe disease is reducing the number of sources or "changing their location" and so on.

OK, so, what we are really saying in this case is that we are recording the activity of a small number of components within the brain with 20 electrodes at the surface of the brain. We enter all these 20 time series in ICA, and ICA returns:

The Independent Components (the elementary waveforms whose linear combinations generate all 20 original recordings)
The "loadings" for each electrode. Or in other words, the mixing levels, how much of each of the elementary components is required to reconstruct the signal at each electrode.

Essentially now we have the signals of the sources and we are re-constructing the signals at the electrodes.

Since we are now in full control, let's think about what happens with just one elementary waveform, let's call it "IC1". If we now plot the loadings on the positions of the electrodes (and interpolate in between) we get a topogaphic map of the strength of the Independent Component across the surface of the brain and to an extent, we can localise its activity. The same can be done for the rest of the components and this is how these topographic maps are calculated.

In the case of EEG and assuming a "standard" montage, it is possible to generate these general images. In the case of Magnetoencephalography (MEG), you can get the actual positions of the electrodes from the manufacturer to generate these plots accurately.

For more information please see this,this and this papers while for a more detailed overview of the method please see this link. The papers are not "de-facto" papers, but give a good overview of the use of ICA in EEG, both for preprocessing but also for analysis.

Hope this helps.


Deriving of magnitude, phase response and impulse response of transfer function

HINT: The frequency response has the form
$$H(e^{j\omega})=a+be^{-j\omega}+ce^{-2j\omega}+be^{-3j\omega}+ae^{-4j\omega}\tag{1}$$
which can be rewritten as
$$H(e^{j\omega})=e^{-2j\omega}\big[ae^{2j\omega}+be^{j\omega}+c+be^{-j\omega}+ae^{-2j\omega}\big]\tag{2}$$
Now note that the term in brackets is purely real-valued. I trust that you can take it from here.


What is the impulse response for 2D functions

As pointed out by Jason R. in his comment, you can easily apply an impulse $\delta(m,n)$ (with $\delta(m,n)=1$ for $m=n=0$, and zero otherwise) to the systems and see what the output is. However, for those two examples the response to an impulse does not say much about the general behavior of the system.

There is an important class of systems called Linear Shift Invariant (LSI) systems which are completely characterized by their impulse responses. Linearity means that if $y_1(m,n)$ is the response to an input $x_1(m,n)$ and $y_2(m,n)$ is the response to input $x_2(m,n)$, then the response to the input $ax_1(m,n)+bx_2(m,n)$ is given by $ay_1(m,n)+by_2(m,n)$. If you check this condition you will find out that your first system is not linear.

Shift invariance means that if $y(m,n)$ is the response to $x(m,n)$ then $y(m-k,n-l)$ is the response to the input $x(m-k,n-l)$. Your second system does not satisfy this condition.

So none of the two systems is LSI and therefore neither of them can be characterized by their impulse response.


Fourier transform of the magnitude of the fourier transformed signal

If it were the squared magnitude, the result would be the auto-correlation function, apart from a time inversion and possibly a scaling factor. If it's just the magnitude, then the result can be characterized as the only even function that has the same Fourier transform magnitude as the original function.


Can I measure non-linearity of a line echo path?

excellent, you have PAM-4!
Send a long sequence of your least-amplitude symbol, and integrate the received power to calculate the energy of your received signal.
Do the same with the three other symbols.
Now, you can solve the equation
$$E_{RX} = \alpha E_{TX} + E_N + F(P_{TX}),$$
with $\alpha$ being the channel gain, which has to be the same for all 4 of your measurements, $E_{RX},E_{TX}$ being the transmitted and received energies (hint: normalize them to the least-power symbol), $E_N$ is noise energy, i.e. noise power times number of samples (which becomes asymptotically constant if you just average enough symbols), and $F$ is the nonlinear error term.
With $E$ known, and $\alpha$ and $E_N$ overdetermined determined from four measurements, $F$ is the reason for the error.