What are Simulated B-mode and phantom images?

B-mode is the most elementary way to form an image from the Ultrasound modality. You basically shift the beam of an A-Mode, either mechanically or electronically. To simulate B-mode, you need a three dimensional model of tissue density and a sound propagation model. There are basically two methods to achieve it, you either do ray-tracing (if the wavelength is much smaller than the distance) or you solve the wave equation directly. In either of these cases you simulate ultrasound propagation through the tissue and re-construct the image that a real ultrasound machine would give you.

A "phantom", more generally, is an item that has been fully characterised and its measurements are known. It is used for calibration. So, for example, you know that Water should return a zero in the Haunsfield scale and Air should return -1000. If you scan a phantom and the image does not agree with its specifications, the machine needs to be re-calibrated / repaired.

The simplest example of a phantom in "everyday life" is that gray-scale ribbon that you see in some screen calibration utilities that direct you to tune your monitor until you can see all the bands in the ribbon, or until some patterns appear / disappear.

A software phantom is exactly the same thing but it is used to "calibrate" or assess an algorithm. You are probably talking about the Shepp-Logan phantom which contains a number of fully characterised ellipses.

If you are working in tomographic image reconstruction then you can compare the image that your technique reconstructs against the "real" image and make inferences.

To create these images, all that you have to do is create an image with shapes of known characteristics (as you can see from the Wikipedia page on the Shepp-Logan phantom).

As far as the "...formula for signal-to-noise ration and contrast-to-noise ratio..." is concenred, it can be found on Wikipedia. It really depends on what you are trying to measure and how and yes, you will need a phantom for that. For more information please see this link.

Image quality assessment is a field on its own in medical imaging.

Hope this helps.


What are the types of noise (FM,AM, MSK modulated) that I can use for jamming GSM signals?

Almost any kind on in-band interference (like Gaussian noise) will mess up the GSM signal and prevent its correct reception. Note, however, that:

Transmitting on the GSM band without a license is illegal in many locations.
Tampering with cell phone operation in general is also illegal.

If you insist on doing this, one approach would be to set up a GSM base station using OpenBTS in an environment where your signals won't interfere with other base stations and won't be otherwise detected (for example in a laboratory with RF shielding). You can use unlocked phones to connect to your OpenBTS base station, and then measure the effect of your interefernce on their performance.


What is the maximum data rate a wireless link can support using BPSK for given BER and SNR?

In BPSK, BER is equal to $Q(\sqrt{2E_b/N_0})$, so you need $E_b/N_0>11.3$ to have a BER no larger than $10^{-6}$. Since the bandwidth is 200 kHz, the maximum rate is 400 kb/s. So, since $SNR\approx136$,

$$\frac{Eb}{N_0}=\frac{S}{N}\frac{W}{R}=136\frac{200\cdot10^3}{400\cdot10^3}=\frac{136}{2}=68$$

Since $68>11.3$, you can be sure that you can transmit at 400 kb/s and meet your BER requirement.

Note that I'm assuming baseband transmission (so that the bit rate is twice the bandwidth). Note also that the question, as posed, does not make a lot of sense, since the bandwidth will essentially determine your bit rate.


Why Successive Approximation, Analog to Digital converter is most commonly used?

For many applications, the technology often selected is the lowest cost one that meets the most common requirements (audio sample rates, etc.)


There is a mysterious factor when compute coeffiencts

It sounds like you're comparing the coefficients for two different filters:

A fourth-order Butterworth filter designed for cutoff frequency $f_c$
A cascade of two second-order Butterworth filters, each designed for cutoff frequency $f_c$.

The two are not equivalent to one another, although in the passband and stopband you will see similar magnitude responses. Specifically, at the cutoff frequency $f_c$, you'll see 6 dB of attenuation from filter #1, while filter #2 would give you 12 dB. Therefore, it's not unexpected that you would get slightly different coefficients for the two implementations.


Independently of what they call it, the process described in the paper is precisely the Fast Fourier Transform. This is made obvious by the FFT algorithm they describe, including the Gold-Rader bit-reversal procedure.

On a side note, if one tries Googling "SFFT", the S in the acronym stands for anything from sliding to spherical, but the F's always stand for Fast Fourier.


Why when you use any filter, it results first a steepness at the beginning to a certain level?

If the offset is known, or can be determined before applying the filter, you can subtract the offset from the signal before applying the filter, then add the offset back after filtering.  If you don't know the offset, neither will the filter when it first starts.


Frequency at which the maximum amplitude occurs of a signal

What you want in this case is to find the argmax of the FFT, the index (function 'argument'/arg) which yields the maximum (max), which is obtained via the findpeaks() function in matlab. If you're working in a different language, just search for argmax algorithm. It's not so easy to find the argmax (although the derivative of the signal will tell you a great deal).

The argmax given via findpeaks() will give you the index into the FFT at which the peak amplitude occurs.

Note: Using this technique will only give you a peak amplitude for a frequency $$\pm \, \frac{Fs}{2 \cdot \text{NFFT}}$$

But first off, you need to take the real part of your signal, just look at the second half of it. It's symmetrical.

When you get the argmax, you can compute the frequency of the peak amplitude by taking the argmax and multiplying by the width of your frequency bins, which is $$\frac{Fs}{\text{NFFT}}$$

For example:

If $$argmax(FFT(x[n]) = argmax(X[n])$$

then $$f_{peak} = \frac{argmax(X[n])*Fs}{\text{NFFT}}$$

where $$Fs$$ is the sampling frequency and $$\text{NFFT}$$ is the size of the FFT.


Question regarding DTFT of a complex signal

I think you almost have the solution. As you suggest in your comment, if you express your signal as :
$$
x(n) = j\cdot \delta(n) - j\cdot \delta(n-4)
$$
As DFT is linear, you will get :
$$
\begin{align}
X(k) &= \sum_{n=0}^{N-1}x(n)\cdot e^{-2j\pi \frac{kn}{N}}\\
 &= j\cdot \sum_{n=0}^{N-1}\delta (n)\cdot e^{-2j\pi \frac{kn}{N}} - j\cdot \sum_{n=0}^{N-1}\delta (n-4)\cdot e^{-2j\pi \frac{kn}{N}}\\
&=j-j\cdot e^{-8j\pi \frac{k}{N}}\\
&=\text{ }...
\end{align}
$$


Examples where non-linear phase filters are used

Two examples:

Audio: The human ear isn't very sensitive to nonlinear phase (probably because the world filters sound with nonlinear phase filters).  It is, however, sensitive to percussive sounds that "pre-ring"*; i.e., that start making sound before the main "bang".

So for high-fi audio, one often wants to use minimum phase filters (which tend to have nonlinear phase response).  In particular, the use of IIR filters is common (or at least commonly advocated).

Control Systems: In a control system, absolute phase shift (or delay) rather than nonlinear phase shift is the enemy.  FIR filters that have make nice pretty-looking output pulses in response to input pulses make systems that cannot be tuned to the highest possible bandwidth without oscillating -- and systems with pure delay (as you get from a FIR filter) are particularly hard to tune.

So with exceptions (i.e., comb filters can be used to advantage, if you're careful) control engineers want any filtering done with IIR filters, and they tend to want that to be kept to a minimum.

* Or at least there are audio equipment designers that believe this -- I'm parroting things I've read.


Seamless audio loops via Fourier transforms?

Try using some sort of pitch estimator to determine the pitch of your sample.  Or hum in tune to a note of a known musical pitch (from a pitch pipe for instance).  Then overlap and cross-fade frames, not 4k apart, but a shorter distance apart that is an integer multiple of the period corresponding to that pitch.


To approximate the volume ($V_{counted}$) you need to count all the voxels.

To roughly approximate the surface area ($A_{approx}$) count all the voxels that have an "empty" voxel as a neighbor. (Have a look at Marching cubes or Marching tetrahedrons for more inspiration and a more detailed discussion on how to determine a surface element).

In order to determine the sphericity I would calculate the minimal bounding box$^1$, take the largest axis and regard it as the diameter: $V_{calc} = \frac{4}{3}\pi r^3$

Since your features ($f_{vol}$, $f_{area}$, $f_{sphere}$, $\dots$) should be in the interval $[0, 1]$ (or 0% - 100%) you simply need to divide the smaller value by the larger value:

$f_{vol}(V_{counted}, V_{ref}) = \frac{\min(V_{counted}, V_{ref})}{\max(V_{counted}, V_{ref})}$

$f_{area}(A_{approx}, A_{ref})$ and $f_{sphere}(V_{counted}, V_{calc})$ accordingly.

$f_{compactness} \tilde{=} f_{sphere}$, since a sphere is the most compact volume.

I hope you can understand, what I am trying to say, and hopefully the answer you are looking for is somewhere in here.

$^1$ if your data is not aligned, you may have to align it first (PCA).


Assuming that the tones are sufficiently sinusoidal, here's a possible scheme, with a bandpass filter being controlled by a frequency estimator:

If the tunable bandpass filter's output is at about the same power as the input signal, then the input consists mainly of a sinusoidal tone. It is perfectly fine if the frequency estimator is a simple one that only gives a meaningful output for the tone to be detected.


Depends on your algorithm. Assuming overlap add/save and a filter length of N we would get:

Zero pad to 2*N and forward FFT: 4*N*log2(2*N) multiplies
Spectrum multiply: 8*N multiplies (4*N if input and filter are real)
Inverse FFT: same as 1
Overhead for padding, overlap handling, house keeping etc.

There are between 1 and 1.5 times the number of adds. There are also various ways to optimize this.

In general, convolution is not always the best way to create reverb. It is computationally expensive, it adds a lot of latency, and it is very difficult to update or modify in real time.


Parks-McClellan - compiling C++ code from Iowa Hills page

You need to back track on that site.
There is a complete set of source code and header files available that you will need.

This is all about the FIR-Filters.

It links to this page with examples.
You need the "Code Kit Download."
This will get you the FIRFilterCode.h.
Main.h will be missing because that belongs to the GUI  - which they don't seem to provide code for.
However, it should be fairly simple to fix.
Remove the include Main.h, then remove all references to TopForm.
You'll probably want to replace them with some kind of logging since all they do is provide feedback to the user.
The involved lines are 273+274 and 278+279.
Notes to that effect are included in the comments in NewParksMcClellan.cpp


Averaging AR-coefficients is a risky business - this could cause drastic shifts in the position of the zeroes of the corresponding polynomial, resulting in instabilities!

If I were to approach the problem the first thing I would try would be the following:

Compute the spectrogram (short-term fourier transform)
Run VQ on the spectral slices with a codebook size of 16, using euclidian distance on the log
Once you have 16 spectral template, convert them back to autocorrelation functions, solve yule-walker for the 16 templates and you get your 16 models.

Or another "safe" way:

For each frame, extract the AR coefficients.
Convert them into reflection coefficients or any other alternative representation like line spectral pairs which can be safely interpolated.
Run VQ on these with a codebook size of 16 to get your 16 models.
Convert back to AR coefficients.

But I really doubt this is a path worth exploring. My two intuitions would be that:

If an audio file is long and diverse enough the 16 specialized models you would get out of it won't be very different form the 16 "universal" models you would obtain out of running the procedure once and for all on a large collection of audio files.
Learning 16 models from a large collection of audio files is likely to give you the kind of naive "fixed predictors" used in FLAC (polynomial interpolators).

One related result is that running a KLT on mel-frequency spectra on a sufficiently long audio clip is going to yield something very similar to the DCT base - that's why audio codecs like mp3 or AAC are not using per-file optimized transforms - the DCT is a good, universal one, and the advantage of using universal transforms is that they don't even require any side-information. This is what motivates my intuition - you are looking at a large dataset (an entire audio file), 16 is small. You'd rather think of 16 good universal predictors and use them than instead of trying to optimize your 16 predictors per file.


What is the minimum number of time points recommended for accurate instantaneous phase estimates using the Hilbert transform

Forgive me for not having a numerical answer to the question. But here are a few things to think about. The answer depends on how the Hilbert transform (HT) is implemented. If the HT is implemented using an $N$-point tapped-delay line filter-like block diagram the first and last $N$-1 output samples are incorrect. So in this situation four input times samples is definitely not enough samples to produce useful results.

If HT is implemented by way of a forward FFT, zero the neg-frequency spectral components, and inverse FFT (as is done in Matlab) you have a better chance of computing meaning results. However, in this second implementation the input time sequence must be long enough to contain at least one cycle of all the spectral components of your input signal in order to compute somewhat meaningful results. It seems to me that four times samples is too short to reliably estimate any signal parameters of most real-world (information-carrying) signals.


camera calibration change orientation of axis

The question is a bit unclear (and am sorry to say, the code even more so). Assuming your 3d-2d point correspondences are correct, all you'd need to do after solving the DLT is to rearrange the orientation of the axes while leaving their directions and the parity of the frame unchanged (i.e. ensure that the determinant of the rotation matrix is +1).


Why could window function help to get more accurate specific frequency amplitude?

As promised, here is a little derivation to give some insight into the effects of windows on DFT values done in reverse.
Suppose that I have a DFT from some signal.
\begin{equation}
\begin{aligned}
X[k] &= \sum\limits_{n=0}^{N-1} x[n] e^{-i\frac{2\pi}{N}kn}
\end{aligned}
\end{equation}
I look at it, decide it needs a little "smoothing", so I decide to average neighboring DFT values to make a new DFT bin.  Let's call it $Y$.  The weights I'm going to use are just a "twist" of the (rescaled) VonHann values.  If you plug in $\omega=0$ this becomes clear.
Plug and chug.
\begin{equation}
\begin{aligned}
Y[k] &= ( -e^{i\omega}, e^{i\omega}+ e^{-i\omega},-e^{-i\omega} ) \cdot (X[k-1],X[k],X[k+1]) \\
&=  -e^{i\omega} X[k-1] + ( e^{i\omega}+ e^{-i\omega} ) X[k] -e^{-i\omega} X[k+1] \\
&= \sum\limits_{n=0}^{N-1} \left[ -e^{i\omega} e^{i\frac{2\pi}{N}n} + e^{i\omega}+ e^{-i\omega} -e^{-i\omega} e^{-i\frac{2\pi}{N}n} \right] x[n] e^{-i\frac{2\pi}{N}kn} \\
&= \sum\limits_{n=0}^{N-1} 2 \left[ \cos ( \omega) - \cos \left( \frac{2\pi}{N}n + \omega \right) \right] x[n] e^{-i\frac{2\pi}{N}kn} \\
\end{aligned}
\end{equation}
Now, my weighted average has all the sudden become a window function.  Again plug in $\omega=0$ and you can see it is the VonHann.
Time for one of those trig identities.
\begin{equation}
\begin{aligned}
\cos(A) - \cos(B) &= \cos\left( \frac{A+B}{2}  + \frac{A-B}{2} \right) - \cos\left( \frac{A+B}{2} - \frac{A-B}{2} \right) \\
 &= -2 \sin\left( \frac{A+B}{2} \right)  \sin\left( \frac{A-B}{2} \right) \\
\end{aligned}
\end{equation}
Substitute it in:
\begin{equation}
\begin{aligned}
Y[k] &= \sum\limits_{n=0}^{N-1} 4 \sin \left( \frac{n}{N}\pi \right) \sin \left( \frac{n}{N}\pi + \frac{\omega}{2} \right)
x[n] e^{-i\frac{2\pi}{N}kn} \\
Y[k] &= \sum\limits_{n=0}^{N-1} w_{\omega}[n] x[n] e^{-i\frac{2\pi}{N}kn} \\
\end{aligned}
\end{equation}
When you plug in $\omega=0$ you see why the VonHann is also called the sine squared window.
So, here is a family of window functions whose base member is the VonHann.
\begin{equation}
w_{\omega}[n] = 4 \sin \left( \frac{n}{N}\pi \right) \sin \left( \frac{n}{N}\pi + \frac{\omega}{2} \right)
\end{equation}
Like I said, I'm not a window expert as I don't care for them much, but I don't think you'll find these on the list in the references.  Trig equations are slippery though, it could be one of those in different form.
Still it is interesting.
Hope this helps.

For those who say it might be worthless, I say not so fast.
Let $\omega = \frac{2\pi}{N}$
\begin{equation}
w[n] = 4 \sin \left( \frac{n}{N}\pi \right) \sin \left( \frac{n+1}{N}\pi \right)
\end{equation}
Which means is it zeroes out the last sample as well as the first and shifts the center of the window to $(N-1)/2$, the center of the sample set, vs $N/2$ which is the center of the repeat frame.

Extra credit:
Which value of $\omega$ cancels the "The twisting caused by being off bin." which is point #2 in the reference in my first answer.


What are the advantages of LTI ( Linear Time Invariant ) systems over other systems?

It's important to realize that in practice many types of systems are used, and only some of them can be regarded as (approximately) LTI. The (didactical) advantage of treating LTI systems in a basic signals and systems course is the elegance and relative simplicity of the underlying theory. Stability and causality are easily checked, and the input-output relation is conveniently described by convolution (in the time domain) or multiplication (in the frequency domain). The Fourier transform is a powerful tool for analyzing LTI systems.

Another advantage is the relative ease with which such systems can be designed. Just think of the vast number of filter design methods.

Even though LTI systems are often used in practice, you will find many situations where other systems are needed. A simple example of a linear but time-varying system is a modulator, which multiplies a signal with a given function. A decision circuit in a receiver is a non-linear system. An adaptive equalizer is a time-varying and non-linear system, since the adaptation of the coefficients depends on the input signal. It would be easy to continue this list of practical examples of non-LTI systems. You will probably learn about such systems in more advanced signal processing courses. However, a deep understanding of LTI system theory will help in grasping these more advanced concepts.


What is issue in square wave reconstruction?

Well this goes to show that Fourier series is just approximation that gets more and more correct when you add more harmonics. Take a look at this:

$\dfrac{4\sin\theta}{\pi}$ is just first harmonic.

Harmonics are integer multiple of base frequency as you can see: $\sin3\theta$, $\sin5\theta$ etc. And the more of them you add to the first harmonic the more your final signal will look like a square wave.

The square wave you showed on the left will have infinitely many harmonics, and the one on the right is just a sum of first few.

GIF author Cmglee


Library/strategy for rotate & scale image alignment/registration

One option is keep all the edge points, or to identify key points (junctions, corner) on both images, and then perform a (semi-)rigid matching, such as Coherent Point Drift, which is robust to missing points.

Further, you could split the biggest images into smaller overlapping patches, and adapt the template to the different patches.


Apply butterworth filter (lowpass) to a signal

I found a way to fix the issue by changing the calculations, so it seems like I had the wrong Nyquist frequency and therefore the wrong cutoff frequency. I hadn't realised that the Nyquist frequency in this case didn't need to apply to the highest frequency present in the sound, but rather to the sampling frequency used to generate the sound (it sounds silly now that I know).

The answer (also in the EDIT section of the question) is:

#Creation of the filter
cutOff = 1000 # Cutoff frequency
nyq = 0.5 * sf
N  = 6    # Filter order
fc = cutOff / nyq # Cutoff frequency normal
b, a = signal.butter(N, fc)

The .filtfilt method didn't work because I had the wrong filter.


You may have a look on Paul Zarchan's lecture note about tracking a Sine Wave - Paul Zarchan - Fundamentals of Kalman Filtering: A Practical Approach - Tracking a Sine Wave.

You also may have a look on Tracking a Sine Wave with Kalman Filter - How to Account for Offset (DC Signal)?

Yet in practice, as @gsmafra wrote, using the DFT, since you know the frequency, is the best choice SNR wise.


Complex band shifting is a little more flexible than real band shifting and filtering, and may be a little more efficient too, though if you change the sample rates when moving a signal to baseband that may not be true.

You are correct that you can move your signal to 0 to $f_2-f_1$ Hz in the method that you describe.  You could do the same thing in one step by multiplying the pass band signal by a complex sinusoid with a frequency of $-f_1 $.

Also, when you band shift with a real sinusoid you have to watch out for aliasing.  Band shifting with a complex sinusoid never causes aliasing (unless you reduce the sample rate at the same time, of course).


Relationship between filters, Fourier magnitude, and frequency

(This explanation is going to be rough around the edges so that I can build your intuition. There are subtleties involved, but this should be good enough to start you off).

Forget about the dB scale for a moment. Let us take a step back:

Imagine that you have some signal that occupies a bandwidth from say, 50 Hz to 100 Hz. Its spectrum is flat across this band. That is, it has the same power at 50 Hz, as it does at 51 Hz, as it does at 52 Hz, etc etc. In this example, we say that it has a constant (flat) power-spectral-density within this band. Let us say that the value of the PSD in this band is 1 Watts / Hz. This means that across 1 Hz of frequency, we will measure 1 Watt in power. Across 50 to 100 Hz, we would measure of course, 50 Watts in power. (To compute your PSD values, they are simply the absolute value of your fourier amplitudes, squared).

So far so good.

Now, imagine a filter, that removes frequencies above 75 Hz and also removes frequencies below 50 Hz. It however keeps frequencies between 50 and 75 Hz. To "keep" a frequency, means it multiplies the power of a signal at that frequency by 1. To "remove" a frequency, means it multiplies the power of a signal at that frequency by 0.

So, you have a signal that is 1 Watt/Hz from 50 to 100 Hz. And you pass it through a filter that multiplies all frequency powers from 50 to 75 by 1, ("keeps them"), and it also multiplies all frequency powers less than 50 and above 75 Hz by 0. ("removes them").

If you pass your signal through this filter, what might you get? Well, you just removed all the energy of your signal from 75 to 100 Hz. Previously those frequencies had a power spectral density of 1 Wattage/Hz, however now they are 0. You have just filtered your signal.

If I asked you at this point, "Plot for me, the PSD of your signal after your filtering", you would draw 1's for frequencies in the 50 to 75 Hz range, and draw 0's for all other frequencies. Now if I said, "I dont like this linear scale. I want you to draw this for me on a dB scale. A log scale." This means that you now simply take the log of the PSD axis. So everywhere we had 0 watts/Hz, they are now -infinity watts/Hz dBW/Hz, (db Watts), and everywhere you had 1 Watt/Hz, we now have 0 dBW /Hz.

Taking it so dB scale doesn't change your frequency axis, only your fourier amplitude/power axis.


I would say, in general, that Dynamic Range is a property of the measuring device while contrast is a property of the signal measures assuming it is withing the Dynamic Range of the measuring device.

Let's have an example.
Imagine our dynamic range is measures as 8 bits system (We can derive that from having a sensor with ~48 [dB] DR).
So we have values in the range $ [0, 255] $.

Let's assume we measured 2 different scenes, one spans values in the range $ [4, 232] $ while the other in the range $ [120, 130] $.

So we have a sensor with Dynamic Range of 48 [dB] and the 1st image has high contrast (Big ratio between the high value and the low) while the 2nd image has low contrast.

For instance, the following code in MATLAB:

mA = repmat(linspace(4, 232, 200) / 255, [200, 1]);
mB = repmat(linspace(120, 130, 200) / 255, [200, 1]);

subplot(1, 2, 1); imshow(mA); title('High Contrast'); subplot(1, 2, 2); imshow(mB); title('Low Contrast');

Yielded this result:

Answering Comments

  the upper limit is determined by saturation and the lower limit by noise"

One nice abstraction to Pixel in CCD / CMOS is a bucket.
The light is like rain and the output value - Amount of rain in the bucket.
So when the bucket is full, saturated, we define it as white. Let's say the bucket can hold 256 drops. So if the bucket is full it means at least 256 drops fell into it )Some might spill over).
The minimum number of Drops which we can measure is 1 drop. So we have 256 / 1 DR -> 8 Bit of DR.

  When an appreciable number of pixels in an image have a high dynamic range, we can expect the image to have high contrast. Conversely, an image with low dynamic range typically has a dull, washed-out gray look

Just as can be seen in the figure from MATLAB, when we have saturates buckets, which are white ans we have almost dry buckets (Black) we have high contrast.
Yet if all buckets are filled with similar amount of rain the image indeed looks dull as can be seen.


Sampled signal reconstruction using matlab

I'm not exactly sure what you mean by

  I want to reconstruct the sampled signal

The term reconstruct has a special meaning in DSP and is related to converting a signal from discrete form to continuous using a DAC and a low-pass filter.

I'm going to guess that you are just looking to plot the signal in time domain, and if that is the case then your code would look like this:

[y, fs, nbits] = wavread('signal.wav')
sample_period = 1/fs;                % sample period is just the inverse of the sample rate
sample_time = size(y)/fs;            % total time represented by the samples
time = 0:sample_period:sample_time;  % create the time array
plot(time,y)


What can be a limit to the maximum possible power of a pilot signal

There is no limit other than what is imposed by the FCC and underlying hardware. Any commercial/consumer device has an operating power defined either on the device or in the manual. The max power of the signal (and pilot) will be within those bounds.


Here is how the algorithm works:

So the steps are:

Input -> apply some window (hann in the frame) -> FFT -> Get amplitude > randomize the phase -> IFFT -> overlap and add -> out

Do it in a loop based in the strech Factor.

Paul wrote a Python code for all steps in this picture.


I think the correct expression for the power is
$$ P = \lim_{T\rightarrow\infty}\frac{e^{-aT}-e^{aT}}{2aT}. $$
When $T \rightarrow \infty$, $e^{aT}$ tends to infinity "faster" than $1/T$ and the limit does not exist.
However, an interesting thing happens when $a=0$. The signal $x(t) = e^{-at}$ becomes $e^0=1$, which is a power signal (with power equal to 1 watt).
I'm sure that (using L'Hopital or some other trick) it can be shown that $P=1$ when $a=0$, but I have not tried it. In any case, I would say that the premise of the question is wrong, and in fact $x(t)$ is a power signal for at least one (real) value of $a$.


How to compute fundamental frequency from a list of overtones?

Nope. Difference between overtones is a good point to start, i,e F3-F2, F2-F1. The differences should be all the same or multiple of each other. The smallest one is often the fundamental. It gets more tricky of the spectrum is "sparse", i.e. a lot of the harmonics are missing. Then you need to find a largest possible divisor that turns all frequencies into integers or, to be precise, so that the ratio of frequency to fundamental is within the measurement accuracy of the nearest integer.


Solve for minimum and maximum value of dct coefficient

The max value for a DCT coefficient with this formula would be when $k_1 = 0$ and $k_2 = 0$.  Using this give us the $\cos(0) = 1$ so you simply sum up the pixel values.  Since the pixel values are $0-255$, then the max value is $64*255 = 16320$.
The min value can be computed by putting a pixel value of $255$ whenever the product of the cosines is negative, and a zero when it is positive.  You would have to iterate through the $k_1$ and $k_2$ values to find the smallest value.


I cannot say much about how can we increase the accuracy of your blob detection by your current method till i exactly understand method and implementation; However, i feel it is quite certain that color (alone) quite often vulnerable to false detection due to the fact that many things in the world corresponds to closer to skin color.

The best best, i believe is to use more information - hopefully an independent feature. For example, you can use shape - and apply the classifier who can detect hand/no-hand based on most possible counter shapes. some times the hand shape may have a round blob - but most often the fingers etc would bring the uniqueness necessary for classification.
you can run across many possible scenarios to create a good model.

Think about it, if you would have given this task to a human being - how would they have guessed it?


Unexpected Result When Using Sinc Interpolation

Since Sinc based Interpolation requires you to know the data at any point it can not be done.

You might do a Truncated Sinc Interpolation.
The artifacts you're seeing can be caused by a kernel which is too short or the parameters aren't good.

In order to create a good Sinc kernel you need to know things about the Band Width of the signal and the Sampling Rate, did you took those into account?


Get sound from matlab and convert with FFT

I use this function:

function audio_recorder = MicFFT(Fs, update_rate)

recObj = audiorecorder(Fs, 16, 2);
get(recObj);

% Record your voice for time specified by update_rate.

while 1
    recordblocking(recObj, update_rate);
    myRecording = getaudiodata(recObj);
    subplot (2, 1, 1);
    plot (myRecording);
    xlabel('Time')
    ylabel('Amplitude')
    title('Time domain representation of y(t)')

    L = length(myRecording);
    NFFT = 2^nextpow2(L); % Next power of 2 from length of myRecording
    Y = fft(myRecording,NFFT)/L;
    f = Fs/2*linspace(0,1,NFFT/2+1);

    subplot (2, 1, 2);
    plot(f,2*abs(Y(1:NFFT/2+1)))
    drawnow;
    axis([0 Fs/2 0 0.01])
    title('Single-Sided Amplitude Spectrum of y(t)')
    xlabel('Frequency (Hz)')
    ylabel('|Y(f)|')

end

end

Fs is the sampling frequency whilst update_rate is the length of each recording


How to prove this theorem about the Z transform and final value theorem?

Let me show you a simple way to see this property. Assume $x[k]$ is a causal sequence and let

$$x[\infty]=\lim_{k\rightarrow\infty}x[k]$$

be finite. Then the sequence $x[k]$ can be written as

$$x[k]=x[\infty]u[k]+y[k]\tag{1}$$

where $u[k]$ is the unit step sequence, and $y[k]$ is a causal sequence that decays to zero as $k\rightarrow\infty$. Taking the $\mathcal{Z}$-transform of (1) gives

$$X(z)=x[\infty]\cdot\frac{z}{z-1}+Y(z)\tag{2}$$

The poles of $X(z)$ are determined by the two terms on the right-hand side of (2). Since $y[k]$ is a causal decaying sequence, its $\mathcal{Z}$-transform $Y(z)$ must have all its poles inside the unit circle. Since the first term on the right-hand side of (2) can only contribute a single pole at $z=1$, $X(z)$ cannot have any poles outside the unit circle. If $x[\infty]\neq 0$, $X(z)$ has exactly one pole at $z=1$ from the first term on the right-hand side of (2). If $x[\infty]=0$, i.e. $x[k]$ decays to zero as $k\rightarrow\infty$, then that first term disappears and $X(z)$ has no pole at $z=1$. This verifies the claim: if $x[k]$ has a finite limit as $k\rightarrow\infty$, $X(z)$ cannot have any poles outside the unit circle, and, if $x[\infty]\neq 0$, it has a single pole at $z=1$ in addition to its poles inside the unit circle.

Also note that the final value theorem directly follows from (2). Multiplying both sides of (2) by $z-1$ gives

$$(z-1)X(z)=x[\infty]z+(z-1)Y(z)\tag{3}$$

Taking the limit $z\rightarrow 1$ yields the final value theorem of the $\mathcal{Z}$-transform:

$$\lim_{z\rightarrow 1}(z-1)X(z)=x[\infty]\tag{4}$$

Note that the limit $\lim_{k\rightarrow\infty}x[k]$ must be guaranteed to exist in order for (4) to make any sense.


Software to register a 2D image to a 3D volume

The software Amira from Visage Imaging has some excellent built-in registration features. I never tried that exact problem, though.

It is unfortunately pretty pricy, but you can download a free trial version of it from the website to give it a try. http://www.vsg3d.com/amira/overview

I also found some well implemented registration algorithms in the software MedInria, an open source and free software specially developed for medical image analysis. http://med.inria.fr


Matlab - Bode plot of Lag Filter + Integrator

I tried it in Octave, there's definitely a glitch like you said. I tried it with "Tustin" instead of "zoh", same result.

However I was puzzled by your high sampling frequency. Your lag controller has a pole at s = -0.01 and a zero at s = -0.1, they are both "slow" while your sampling frequency is really high.

How to fix it :
Change your sampling frequency to something much lower than 12 MHz, like 10 kHz for instance and replace "zoh" by "tustin". "Zoh" should only be used when converting a "process" transfer function to the discrete domain in order to take in account the effect of the DAC, we don't use it for compensators.

The root cause however is probably some quantization issue in the c2d function.

FS = 10e3;
T_sample = 1/FS;

accu = tf([0 1], [1 0]);
lag = tf([10 1], [100 1]);

accu_d = c2d(accu, T_sample, "tustin");
lag_d = c2d(lag, T_sample, "tustin");

figure
bode(lag*accu, accu_d*lag_d);

Edit 1 :

The issue seems to be with the lag controller, not the accumulator. You should try to move the pole and the zero to the left, this should fix your problem since you must keep your sampling frequency really high:

try this (same DC gain, but pole and zero further to the left) :

lag = tf([1 1000], [10 1000]);


Estimating SNR of an Environment where obtaining clean speech is not possible

I think it all depends on the SNR region that you are dealing with. Using a voice activity detector (VAD) seems to be worth an attempt, but consider that VADs fail if the SNR is too low. The most simple approach might be using quantiles of the power distribution. Some lower quantiles will probably not contain power that is caused by the speech signal part while high quantiles will contain power resulting from the speech signal. Choosing the right quantiles might be a challenge, but due to the straight forward approach it might by worth a try...

Furthermore, by reinterpreting your problem to noise floor estimation, you might succeed by using some noise floor estimation algorithm (e.g. Minimum-Statistics by Rainer Martin [1, 2]) and calculating the SNR from that.

[1] Rainer Martin, “Spectral Subtraction Based on Minimum Statistics,” presented at the EUSIPCO, 1994.
[2] Rainer Martin, “Noise Power Spectral Density Estimation Based on Optimal Smoothing and Minimum Statistics,” IEEE Transactions on Speech and Audio Processing, vol. 9, no. 5, pp. 504–512, Jul. 2001.


How to get discrete values from a nosiy data set that have a trend?

From the comments you are most interested in counting the number of clear up areas over a certain duration of time, and you point out that you will have occasional erroneous negative spikes even during the up intervals.

You obviously want to filter but you also suggested little knowledge in this space so I will keep it very simple in what should provide a reasonable result.

My suggestion is to do a moving average over the minimum duration you ever expect for an up interval. Pass your waveform through the moving average, and then on the filtered output set a threshold detector with hysteresis (two thresholds) such that you increment your Up counter when passing the upper threshold and then allow your counter to trigger again once your lower threshold is set. Adjust your threshold levels for best results based on your actual signal cases.


Can this be considered wide sense stationary?

for 1, a random process is WSS if its autocovariance and mean ensemble average don't vary as a function of specific time instance, just lag.  As mentioned before, it's hard to conclude this from simply one realization of the process about the ensemble average.

However for 2, you can at least estimate what the autocovariance and mean are if the process is ergodic (which, unfortunately, you can't conclude with just a single realization).  If the process is indeed ergodic, the ensemble average is equal to the time average in the limit as your interval of average approaches infinity, which is kind of nice with regards to computability.


Detection complex vectors by using correlation

The reason we often use absolute value is that the received signal usually has unknown phase.

Consider the case in which $\boldsymbol{y=x_i}\exp(j\theta)$, where $\theta$ is a phase rotation.  The phase rotation is due not to "fading" but due to the fact that the received phase in communication systems is generally arbitrary due to the fact that the receiver's oscillator does not have the same phase as that of the received signal. (Note I ignored noise for the sake of simplicity, but this doesn't change the basic point.)

Now $\boldsymbol{y^Hx_i} = \boldsymbol{x_i^H} \exp(-j\theta) \boldsymbol{x_i} = \|\boldsymbol{x_i}\|^2 \exp(-j\theta)$.  Depending on the value of $\theta$, the real component of this term may be large or small.  That is what Stanley Pawlukiewicz was pointing out.  In his example, $\theta = \pi$, so the output of the correlation would be $-\|\boldsymbol{x}\|^2$.  The real component would be small, and by using the real component, one would erroneously conclude that the two signals were not aligned, when it fact they were aligned.  Using the absolute value corrects that problem.


What sonic artifacts are caused by phase vocoders?

The term "phase vocoder" may encompass more than one resynthesis method.  One method uses a pruned frequency reassignment of the DFT basis vectors, and thus resynthesizes sound using a sparser set of basis vectors with slightly different frequencies than were present in the original DFT analysis windows.

This sparse basis may not be able to recreate a full spectrum wide-band signal, such as produced by a transient, as well as the original non-reassigned full DFT basis.

Phase vocoder analysis also makes the assumption that the spectrum is fairly stationary between FFT frames, an assumption that can also be violated by transients.

One term for the artifacts produced by this type of analysis/resynthesis is transient smearing.


IIR filter design in digital domain using the magnitude squared

To solve the case that you mentioned...

You have 2 variables to determine, so you need two relationships to resolve the two variables.  I'm going to use $k$ and $a$ as the variables to make this easy to type up.

$$H(z) = k\frac{1 + z^-1}{1-az^-1}$$

Start by considering the passband gain.  Use $f = 0$ for this.  Assume you want unity gain at $f_0=0$.

Assume:  $H(f_0) = 1, f_0 = 0$

Substitute $e^{i2\pi f/f_s}$ for $z$, $f_s$ is your sampling rate, set $f = 0$ and solve for $k$ to satisfy $H\left(f_0\right) = 1$.

From this you get $k = \frac{(1-a)}{2}$

Now work on the gain squared at your desired corner frequency ($f_c$) to determine $a$.

$H(f_c) = -3\textrm{ dB}$ (magnitude squared will be $-6\textrm{ dB}$ as you've stated)

We'll work with the magnitude squared at $f_c$ and set the gain to $1/2$ ($-6\textrm{ dB}$).

$\lvert H\left(f_c\right)\rvert^2 = \frac{1}{2}$

This time substitute $e^{i2\pi fc/f_s}$ for $z$.

To simplify the arithmatic you can solve this equation:

$$
\left(\frac{\lvert H(f_0)\rvert}{\lvert H(f_c)\rvert}\right)^2 = 2
$$

This eliminates the factor $k$.

You will end up with a quadradic relationship in $a$.  Solving for $a$ yields:

$$
a = \frac {1 - \sqrt{1-\cos^2\left(2\pi\frac{f_c}{f_s}\right)}} {\cos\left(2\pi\frac{f_c}{f_s}\right) }= \frac {1 - \sin\left(2\pi\frac{f_c}{f_s}\right)} {\cos\left(2\pi\frac{f_c}{f_s}\right) }
$$


Algorithm to detect if a spoken phrase matches a previously spoken phrase

MFCC is an established way of handling speech recognition. You will find plenty of libraries and tools (e.g. yaafe, SoX) that can extract them for you. Implementing them yourself is also quite easy. Yaafe provides plethora of other audio features you might like to play with.

Mind that MFCC are not very robust against noise, so you might like to apply some noise reduction beforehand if required. Some people advocate raising the log-mel-amplitudes to a power 2 or 3 before applying discreet cosine transform.

Some libraries, like scikit learn, next to all machine learning algorithms provide also DTW. On your place I would just try MFCC and see how it goes from there.


Understanding Eigenface Face Recognition Algorithm Notations

G1, G2... Gm are column vectors of length 65536.  The matrix A has columns that are the principal components.


How to know if a audio file is real lossless using its spectrogram?

Eventhough there could be opposing theoretical answers to your question, practically speaking it may have a solution which depends on the lossy codec you are using. In that case what you eventually need is a practical knowledge of codec under concern that is used for the audio clip. However, considering the number of different audio codecs and the fact that there could be any number of nonstandard implementations to exist, your problem may not have a guarantied solution. It is because for every different lossy auido codec a different procedure (provided in theory it does exist) might be necessary to test whether a given piece of audio belongs to its output or not. And since for a nonstandard/proprietary/unknown audio codec, you will not have such a procedure at hand you may not even be able to test to see it.

For example, mp3 (MPEG-I layer 3) standard defines an optional low-pass filter about 16khz, before encoding the raw pcm audio, hence if you find some signal acitivity above 16 khz inside an audio clip, then it is either not an mp3 file or the mp3 profile used to encode that clip, did not activate the lowpass prefilter (as it can be one of the many optional features of the codec) Note that in either case you cannot say that it is a raw file just because it's not an mp3 as it might come from another lossy codec... That is the problem...

Finally, if a lossy codec will include a distinctive feature as a mandatory (nonoptional) requirement of the standard, then you can check for that property to at least decide whether the given piece may belong to its output or not.


A system that perfoms Fourier Transform operation - is it an LTI system?

The Fourier transform operator $\mathscr{F}$ is a linear one; i.e.,
$$\mathscr{F}\{x(t)\}=X(f) ~,~ \mathscr{F}\{y(t)\}=Y(f)  \implies \mathscr{F}\{\alpha x(t) + \beta y(t) \} = \alpha X(f) + \beta Y(f)$$ And therefore the system that implements it will be linear.

However time invariance, which is tested by the method

$$\mathcal{T}\{x(t)\}=y(t) \implies \mathcal{T}\{x(t-d)\}=y(t-d)$$

does not apply to the Fourier Transform operator, because the input and output functions do not belong to the same domain (namely time here)

However, the physical system which produces Fourier transform of the given signal $x(t)$ can be associated with a time-dependent Fourier transform which computes the function

$$X(\omega,t)= \mathcal{TF} \{ x(t) \} = \int_{-\infty}^{\infty} x(t+\tau) w(\tau) e^{-j \omega \tau} d\tau $$

Then it can be shown that this operator is time invariant in the sense that,

$$ \mathcal{TF} \{ x(t) \} = X(\omega,t) \implies \mathcal{TF} \{ x(t-d) \} = X(\omega,t-d)$$

Lets elaborate on this. Now lets define STFT as
$$X(\omega,t) = STFT \{ x(t) \} = \int_{-\infty}^{\infty} x(t+\tau) w(\tau) e^{- \omega \tau} d\tau $$ It's easy to show that $$X_2(\omega,t) = STFT \{ x(t-d) \} = \int_{-\infty}^{\infty} x(t-d+\tau) w(\tau) e^{- \omega \tau} d\tau = X(\omega,t-d) $$ Hence it's a shif-invariant operator.

EDIT: Below shown is a trivial redefinition of the above STFT, which does not confirm to the standard notation, hence may not apply. I don't delete it for now however. Shift-invariance is shown for the above definition only.

Now a similar form discussed in the provided link is obtained by a change of variable $t+\tau = \tau'$ on the STFT definition integral:

$$ X(\omega,t) =  \int_{-\infty}^{\infty} x(\tau') w(\tau'-t) e^{-j \omega (\tau'-t)} d \tau' $$
and replacing $\tau'$ with $\tau$ yields the equivalent form of STFT as
$$X(\omega,t) = STFT \{x(t)\}= e^{j \omega t} \int_{-\infty}^{\infty} x(\tau) w(\tau-t) e^{-j \omega \tau'} d \tau $$

Now lets show its time invariance by the input $x_2(t) = x(t-d)$ :
$$ X_2(\omega,t) = STFT \{ x_2(t) \} = STFT \{ x(t-d) \} = e^{j \omega t} \int x(\tau-d) w(\tau-t) e^{-j \omega \tau} d\tau $$ Make a change of variabes inside the integral operator as $\tau-d = \tau'$ yields

$$ X_2(\omega,t) = STFT \{ x(t-d) \} = e^{j \omega t} \int x(\tau') w(\tau'-(t-d)) e^{-j \omega (\tau'+d)} d\tau' $$ which simplifies by replacing $\tau'$ with $\tau$ and extracting $e^{-j\omega d}$ outside of the integral

$$ X_2(\omega,t) = STFT \{ x(t-d) \} = e^{j \omega t} e^{-j \omega d} \int x(\tau) w(\tau-(t-d)) e^{-j \omega \tau} d\tau $$

$$ X_2(\omega,t) = STFT \{ x(t-d) \} = e^{j \omega (t-d)}  \int x(\tau) w(\tau-(t-d)) e^{-j \omega \tau} d\tau = X(\omega,t-d)$$

Which proves that this form is also shift-invariant on the $t$ variable.


Cross Power Spectral Density of Unevenly Sampled Data

One way would be to interpolate the signals to produce evenly-spaced samples and then calculate the CPSD as normal.


This question is too broad and without some more information about the application and the constraints impossible to answer. For actual filter "design" a dedicated DSP is probably not a great choice and you are better off with a high end general purpose processor such as an Intel Xeon or then new-core I7s, Xeons also work well in multiprocessor setups.

Filter implementation really depends on whether you have high-complexity/low-clock or low-complexity/high-clock rate. Best solution for an 100k+ tap room response filter for audio, or a line EQ for a cable modem are totally different. In general for low complexity applications ASICs are faster than FPGAs are faster than general purpose DSP processors.


Loop bandwidth for symbol timing recovery

There is an optimum loop BW that maximizes the SNR, and this applies to both Symbol Timing Recovery as you inquire about as well as Carrier Recovery. The specific answer depends on the characteristics of the noise source involved in your system, such as clock jitter and phase noise, and the modulation you are using; specifically how the signal energy is distributed in frequency. Without knowing those an exact answer cannot be determined, but I can at least share some insights from my own experience with it.

To state generally, I have found that the optimum loop BW to usually be somewhere between R/100 and R/20, where R is the symbol clock rate being recovered.  I start with a value in this range and then through simulation with a valid noise and waveform model (or lab measurement) determine the actual optimum using SNR metrics vs Loop BW setting.

As for damping factor, I also start with a damping factor of 0.5 and then through evaluation using a step or impulse response test adjust for possible ringing due to parasitic delays that may be present in the loop.

In addition to the SNR considerations given above, there are also acquisition time constraints or possible dynamic noise conditions that would also motivate you to have as fast of a timing loop as possible: The faster the loop, the faster the acquisition time, and recovery time from an interference event, or the ability to track a highly dynamic change if it exists.  There may also be waveform specific constraints such as the duration of a header in a packet based waveform used to establish timing. However as the BW for the tracking loops increase beyond a significant fraction of the timing symbol clock rate (such as >R/20), tracking errors will start to dominate the noise in SNR, since the loop will start to track out the modulation of interest.

To explain further from the perspective of a trade of Loop BW and SNR, the Loop BW sets the boundary at which your recovery loop will track timing noise (jitter/phase noise) in your signal.  Timing noise has components at all frequencies (such as that viewed with a phase noise spectral density), and the lower frequency components of the timing noise (the change in zero crossing locations that is moving very slowly) will be tracked by the timing loop, and therefore suppressed. However the higher frequency components will not be tracked by the timing loop and therefore remain and contribute to noise as part of SNR. The timing loop can be viewed as a high pass filter on the jitter/phase noise components, passing the higher frequency components (including the modulation signals) and supressing the lower frequencies at a rate depending on the order of the loop. Thus you see the motivation to make the timing loop as fast as possible in the interest of tracking out as much timing jitter as possible. However, some of the signal energy of interest will be part of the "noise" being supressed by the timing loop; and this can be quantified by understanding the spectral characteristics of the modulation used. As the loop BW increases, too much signal energy is also removed, and thus it is clear how there can be an optimum setting for the Loop BW that maximizes SNR.

This hopefully provides insight into the approach and what to look for in the development of the system. In conclusion, I set the Loop BW as fast as possible (typically in the range of R/100 to R/20 where R is the symbol rate) while monitoring through simulation with accurate noise and waveform models the SNR metric; and either choose the point of best SNR, or if other factors such as acquisition time or more stringent dynamic conditions are involved, maximize BW without going below a target SNR metric.


I have a mix signals (2 signals with different dc offset). How to recover another dc value which has been buried by the stronger dc value?

In general, to be able to separate two signals, you need two different mixtures of there signals. It's very similar to solving a system of equations with two unknowns. You need two equations if you have two unknowns. If you simply have a signal which is a sum of two other signals, you cannot retrieve the individual components.


Calculating a homography based on detected lines

I will explain two approaches for this:

1)
One approach would require a line matching algorithm. After matching the lines, you could simply use the end points of the lines in order to compute homography. To achieve that EDLine or LSD based descriptors are recently proposed in OpenCV. Also, hashing and fast matching of them are also implemented. Check out the videos here:

http://www.youtube.com/watch?v=MqMjvSkM39k

http://www.youtube.com/watch?v=naSWTlbg3To

The recent opencv_contrib repository contains the source code to these methods.

In the case that the line end points are noisy, you could then directly utilize the lines to compute the homographies. Such papers would then read:

  Internal Report: 2005-V04 Computing Homographies from Three Lines or
  Points in an Image Pair
  G. Lopez-Nicolas, J.J. Guerrero, O.A. Pellejero, C. Sagues

  Internal Report: 2003-V01 Robust line matching and estimate of homographies simultaneously G. Lopez-Nicolas

  Probabilistic Matching of Lines for Their Homography Taemin Kim, Jihwan Woo, and In So Kweon

2) There is one method specific to the fields are given here:

  "Using line and ellipse features for rectification of broadcast hockey video.", Gupta, Ankur, James J. Little, and Robert J. Woodham
  Computer and Robot Vision (CRV), 2011 Canadian Conference on. IEEE,
  2011.

and

  "Combining line and point correspondences for homography estimation.", Dubrofsky, Elan, and Robert J. Woodham. International Symposium on Visual Computing. Springer Berlin Heidelberg, 2008.

The idea is as follows: Any line, parametrized by its coefficients $\mathbf{l}_i=(u,v,1)^T$ maps to $\mathbf{l}_i^{'}=(x,y,1)^T$ in the other image using:

$$
\mathbf{l}_i^{'} = \mathbf{H}^T \mathbf{l}_i
$$

In this form the equation can be directly plugged into the DLT method:

$$
A_i =
\begin{bmatrix}
-u & 0 & ux & -v & 0 & vx & -1 & 0 & x  \\
0 & -u & uy & 0 & -v & vy & 0 & -1 & y
\end{bmatrix}
$$

The only difference is the normalization, which you will find in the references above.

Adding Ellipses:
Any point $x$ lies on the conic section $C$ if $x^TCx=0$. This gives rise to the transformation relation:

$$
C' = H^{-T}CH^{-1}
$$

The references above also explain how to insert this constraint to the DLT algorithm.

Using ellipses and lines, it is possible to derive a robust projective relationship.


https://ccrma.stanford.edu/~jos/filters/ and specifically
https://ccrma.stanford.edu/~jos/fp/Difference_Equation_I.html
http://www.dspguide.com/pdfbook.htm
https://www.amazon.com/Digital-Signal-Processing-Alan-Oppenheim/dp/0132146355
    (need to buy the book)


If two signals have different distributions that you know and can model, can you use that to separate them?

Statistically, yes. Does that mean you can be certain at any point in time? No. But it allows you to transport data.

For example, if I have a mixed signal composed of signal A added to signal B, and I know that the histogram of signal A is non-Gaussian in the time domain, and that the histogram of signal B is Gaussian in the time domain, is there a way to isolate either A or B using this information?

Not to separate them, but to use both to transport data independetly: For example, you can switch the Gaussian source on and off, and detect the presence of it.

I understand that I can use MLE or some similar method to fit a Gaussian to signal B

Uff, you don't "fit" a Gaussian to something that has a distribution like a Gaussian. I know this sounds like nitpicking, but in these contexts, it's very important to make a distinction between the shape of the signal and the signal's amplitude distribution.

My initial idea was to create a Gaussian filter that matches the distribution of signal B, and use that to separate B from A+B, but I'm not sure how I would implement that in practice, or if that even makes sense.

Doesn't make sense, far as I can tell.
I'm afraid you'll have to look into information theory a bit, and there into what mutual information is, how we model that, and noise.
You can then make general statements on the ability of A to transport data in the presence of B, considering the latter "noise". It also works the other way around!
This directly leads to the successive cancellation decoder, where you first (try to) decide the better receivable variable, then subtract the decision from the observation, then decode the other. For example:
Your A is discretely uniformly distributed on ${-2, -1.9, -1.8, \ldots , -0.1, 0, +0.1,\ldots +2\$, your B is Gaussian with zero mean and variance 0.01. You observe  Y = A+B = -1.52; you decide that by far the most likely value for A is -1.5, and with that knowledge you can infer that B must have been -0.02.
You needed distribution knowledge of A and B here to decide which of the two to decode first, and then to make a maximum-likelihood decision on A.
This successive method gives you a lower bound for your ability to get data from A+B given their distributions. However, that's not the same as a joint decoding. General, for a maximum likelihood estimator, you'd need to know the joint distribution of A and B under A+B; then you can, for any A+B, tell the most likely A and B. There's to my knowledge no general formula for this decoding – it really depends on the distributions. Also note that if A and B are statistically independent, this joint density will fall apart into the product of the individual densities. Note, however, also that usually, both A and B will employ channel coding, so that these likelihoods are both time-variant and depend on other decisions – welcome to the wonderful world of coding theory.


Picking out a signal that appears to be noise inside a large signal

Adaptive filters could be used in situations where you have a reference signal, what adaptive filters basically do is to try modelling the signal that is present in captured signal using reference signal. So, we can subtract and/or extract the required signal.

In the situation mentioned in question, you can use adaptive filter if you have reference for either the signal you are interested in or the noise you want to remove. If you have the noise signal as reference adaptive filters will remove the required signal(signal of interest) from input and only noise will be remaining, now this noise can be subtracted from original input signal to get the required signal(signal of interest).

Same could be done if signal of interest is present as a reference, but in this case the adaptive filters output will be the required output.


Orthogonal basis of signal space and the projection of white noise

Leave time dimension aside for a moment and deal with just one symbol at a time. In a PAM modulation the generated functions are not ortho and not normal. Actually, in a 4-PAM modulation you have four basis pulses: $p_0(t)=a_0 \cdot p(t)$, $p_1(t)=a_1 \cdot p(t)$, $p_2(t)=a_2 \cdot p(t)$ and $p_3(t)=a_3 \cdot p(t)$. The pulse $p(t)$ is explicitly defined just because it makes the math easier and coherent: it is a function that generates a 1-dimensional Hilbert subspace.

The task is about mapping data symbols to signals in the Hilbert space, in such a way that at the receiver side the mapping can be reverted, thus, such a mapping must be bijective. On the receiver side, we must compute some kind of metric over the symbol period to come back to the discrete domain. For a Hilbert space, such metric is the inner product. So, in both cases (FSK and PAM) you end up performing, one way or another, the same computation:

\begin{equation}
a = \int_T f(t) g(t) dt
\end{equation}

And this inner product is performed against each orthogonal basis function (one or more). By choosing more than one orthogonal basis functions, you get more parallel channels (more dimensions) to transmit onto if you wish. Then, choose how to encode the data:

BFSK encodes data by signaling one of those channels or the other. So for one of them, $a \neq 0$ (gets mapped to symbol $s_1$) and for the other (assuming orthogonality) $a = 0$ (mapped to symbol $s_0$). If both channels get signaled, you cannot reverse the mapping.
2-PAM encodes data by signaling, within the same channel, different levels. So that if $a \gt threshold$, symbol $s_1$ is decided. If two levels are sent, you cannot reverse the mapping.
M-QAM uses both methods: different levels and different channels within each symbol period.

In short, data can be encoded in the value of $a$ within the same channel (same dimension, like a scaled basis function) so you further have to define decision boundaries for mapping $a$ to discrete symbols, or can be encoded by which channel (or dimension the orthogonal basis function generates) you signal. This resembles vector-space math because Hilbert space is the space of continuous vectors.

So, for the questions:

Q1. Yes. Harmonic sines within the fundamental period are orthogonal under the inner product. Orthonormal just means the norm for each basis function equals 1.

Q2. No. When it is said that noise is uncorrelated it refers to the fact that AWGN has no memory (time dimension), the noise is already uncorrelated before projection onto any basis. Within a symbol period the noise is correlated, its value is its power over that period.

Q3. When the constellation is said to have M points, it means M points get mapped to M different functions in the chosen Hilbert sub-space. The way you choose those functions to be is not defined. You can choose a M-dimensional space (so you only need to find the orthogonal basis function which fits better) or a N-dimensional space with $M \gt N$ and define different regions in that space for symbol mapping. Also, in QAM, $p_r(t) = cos(\omega t) \ne sin(\omega t) = p_i(t)$, and the inner product is zero within one period.

I avoided to talk about the time dimension (when possible) because IMHO it would not help. Unless you are trying to exploit temporal correlation (viterbi, equalizers, noise shaping, etc...) it does not help and only introduces noise in the explanation. Usually, this kind of analysis is performed on a symbol based period. However, according to Fourier analysis, every function in the Hilbert space can be expressed in terms of amplitude, frequency and phase, and a shift in time is a shift in phase, thus the same logic applies. In real implementations you do not find neither infinite signal levels ($a$'s), since power is constrained, nor frequencies (basis function's shape) which are constrained by bandwidth. Time shift (phase) is not infinite neither since you need to start and stop sometime.


If you are sampling 50 samples/sec your nyquist frequency is only 25Hz, so your test signal is being aliased. You should expect different results if your sampling rate is less than double the highest frequency in your signal.


How to determine where to sample for demodulation of BPSK signals?

I have done something similar to this in MATLAB. In my cause, I used an Early/Late Gate clock recovery method to get an estimate of the offset between transmit and receive symbol timing. This method uses 3 samples per symbol - one at the optimal sample time, one that is 1 sample delayed and 1 that is one sample advanced. This works well for on/off keying signals, since there is no negative portion of the signal, but the methods below I think will work better for bipolar signals.

Another method is the Muller and Mueller Algorithm, which uses 1 sample per symbol. This uses the equation $e_n = (\widehat{y_n}y_{n-1}) - (y_n\widehat{y}_{n-1})$, where $e_n$ is the error , $y_{n-1}$ and $y_{n}$ are the previous and current symbol samples, and $\widehat{y_n}$ and $\widehat{y}_{n-1}$ are the decision outputs (-1 or +1 in your case). M&M requires that the carrier is recovered first.

A third algorithm, which is very similar to the one above, is the Gardner algorithm. In this algorithm, the error signal is $e_n = (y_n - y_{n-T})y_{n-T/2}$, where $y_{n-T}$ is the sample one symbol in the past, and $y_{n-T/2}$ is half a symbol in the past. This requires 2 samples per symbol, but is not sensitive to the carrier.

The M&M detector is looking at the peaks, while the Gardner method is looking at the peak and zero crossing to get timing estimates.

In all cases, you'll take the error estimate, feed it into a loop filter whose output is then used to advance or retard the symbol clock for the next timing estimate. The loop filter can be thought of as either a low-pass filter, a smoothing filter or simply as weighting older estimates less than newer ones. The loop filter controls a) how fast the timing converges and b) how much noise affects the estimates.

Another option is to use an equalizer which will then remove timing delay from the signal. This can also help for very noisy signals I believe.

I have used these resources in forming this answer:
http://www.comlab.hut.fi/opetus/333/2004_2005_slides/CarrierTimingRecovery.pdf
http://mobiledevdesign.com/tutorials/radio_matched_filtering_timing/

This article is especially helpful: http://rfdesign.com/images/archive/0901Litwin32.pdf
Here is a link to Gardner's original paper from 1986: http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1096561&tag=1

Note In general this area is called "clock and data recovery" or "CDR". Specifically we're dealing with "symbol synchronization". These terms should aid you in your search.


HINT:

The only thing you need to know is that any sequence convolved with a (shifted) delta impulse is a shifted version of itself:

$$f[n]*\delta[n-k]=f[n-k]$$

Apply this to the expression $u[n]*(\delta[n]-\delta[n-1])$, and realize that what you're left with is a very simple sequence. (The result of (a) will be helpful.)

Then convolve this very simple sequence with $x[n]$ to get the final result.


What fps can be considered as a hardline for real-time performance? Is there any academic paper that describes this?

There are as many ad hoc "definitions" of "real-time" are there are people or places you hear or read about what it supposedly means. Some people claim there is a "hardline" such as 10uS response time, but I believe there is no academic paper that makes that mistake.

Many people implicitly have an informal mental model that considers a system as being "real-time"
•   if, or to the extent that,
•   its actions are manifest to action observers with latencies that can be related (usually in an application-specific way) to the actions’ perceived currencies—
•   i.e., in a time frame that those actions’ latencies and predictability of latencies have acceptably satisfactory value to their observers (e.g., system users). Think of high frequency trading.

The magnitudes of the time frames are application- and situation-specific—e.g., microseconds to megaseconds—not related to if, or to the extent, that a system is a real-time one.

So a system "operates in real-time" if its actions have low enough latency to satisfy your needs. If 15 FPS is a satisfactory enough rate, then the system is real-time as far as your needs are concerned; otherwise it is not--that's all that matters to you.


If you are only decimating by four, then a CIC filter is not the way to go.  CIC filters have two advantages and one big disadvantage.  The advantages are that they don't require any multipliers, just (large) adders, and they can efficiently decimate by really large factors.  The no-multipliers advantage was a big deal years ago when multipliers were expensive and not always available, but they are cheap now.  And since you are only decimating by a factor of 4, the ability to decimate by large factors does not help you either.

The one big disadvantage of CIC filters is that their frequency response is not flat in the pass-band, so designers usually put a FIR filter after the CIC filter whose primary purpose is to compensate for the CIC filter's non-flatness.

My advice would be to continue decimating with one or two FIR filters, though I would recommend using a better filter than [1 5 10 10 5 1].  Try using Matlab's fdatool or Matlab/Octave's "fir1" command to design a better filter.


Minimum Sampling Rate of Bandpass signal

The problem is that the general criterion for bandpass sampling (which results in $25\text{ kHz}\leq f_s\leq 30\text{ kHz}$) applies for general signals, whose spectrum may not be symmetric around $f_c$. In your example, sampling at 10 kHz would fail because spectrums of different shapes overlap.

The most complete resource I know on bandpass sampling is "The theory of bandpass sampling", by Vaughan et al, IEEE Trans. on Signal Processing, 1991, Vol. 39, Issue 9 (http://dx.doi.org/10.1109%2F78.134430). At the start of page 4, it specifies that in the special case of DSB signals, the minimum sampling interval is $1/B$, where $B$ is the bandwidth (for both sidebands). This corresponds to your case, where the minimum sampling rate is 10 kHz, equal to the signal's bandwidth.


This is a hard problem because the image data is encoded using variable length codes (VLC's).  This means that the starting bit of each DCT coefficient is dependent on all of the previous coefficients.  A single bit error will often produce a different code length for that coefficient and in turn throw off the alignment for the rest of the image.

I have written automated programs before to toggle bits around where the first corrupt block is and try to find where the errors might be.  Basically, it just toggles a bit, redecodes the image, and runs some kind of metric to determine if it helped or not.  This method may or may not work if there are multiple bit errors close together.  The key is that you don't necessarily need to fix all bit errors, just produce a VLC that is the right length to recover alignment.


What is the significance of the linear phase in the phase response of an M-point moving average filter?

The significance is a constant time delay for all frequencies. Time delay is the derivative of phase with respect to frequency, so given a linear phase as shown, the time delay is constant. Notice that the abrupt steps in phase actually only occur when the magnitude goes through zero so does not effect the nature of it being a "linear phase" filter.
An M-point moving average filter will have a time-delay (for all frequencies) of $M/2$, so if you increase the number of points, the phase slope versus frequency will increase.
Why do we like linear phase filters so much? Please see this post on that:
Why is a linear phase important?
Any finite impulse response (FIR) filter with symmetric or anti-symmetric coefficients (The OP's case with coefficients all 1 is symmetric) will be a linear phase filter.


Your goal should be achievable. The key driver for the filter's complexity is the size of the transition band, which you haven't specified. If you can live with a slow rolloff, then you can use a relatively low order filter.

If, as is more likely, you want to isolate the band from 5 to 50 Hz with relatively sharp cutoffs, then as you have started to discover, you will want to use a multi-stage approach. Instead of decimating by a factor of several hundred all at once, you decimate by a series small factors (something like 2, 4, or 8). By gradually decreasing the sample rate, you end up with filter designs that are much less complex.

This previous question has a very similar aim to yours, only with slightly different filter specs and sample rate. You should be able to adapt the information there to your application.


Considering @Stanley Pawlukiewicz's hint, you should take into account that $e^{~j \pi} = -1$, thus $(-1)^n = e^{~j \pi n}$, for any integer $n$.

Therefore, time-domain multiplication by $(-1)^n$ corresponds to a frequency-domain shift by $\pi$ radians in digital frequency.

Recall the Fourier transform's theorem: $\mathscr{F}\{x[n]~e^{j \omega_0 n}\} = X(\omega - \omega_0)$ where $\mathscr{F}\{\}$ is, in this case, the Discrete-Time Fourier Transform (DTFT).

Now, since

$$\mathscr{F}\left\{\frac{\sin (\pi n/2)}{\pi n}\right\} = \sum\limits_{k=-\infty}^{+\infty} \Pi\left(\frac{\omega + 2k\pi}{\pi}\right)$$

then, answering to your question:

$$\mathscr{F}\left\{(-1)^n \cdot \frac{\sin (\pi n/2)}{\pi n}\right\} = \sum\limits_{k=-\infty}^{+\infty}\Pi\left(\frac{\omega - \pi + 2k\pi}{\pi}\right)$$

which is a box function of amplitude $1$ and width $\pi$, centered on $\omega = \pi$.  The box function or rect function is

$$\Pi(u) \triangleq \begin{cases}
0 & \mbox{if } |u| > \frac{1}{2} \\
\frac{1}{2} & \mbox{if } |u| = \frac{1}{2} \\
1 & \mbox{if } |u| < \frac{1}{2}. \\
\end{cases}$$

Since the DTFT is periodic in $\omega$ with period $2\pi$, there is a box centered on $\omega = (2k-1)\pi$, for every integer $k$.


Determine modulation function from spectrum?

The Fourier transform is reversible, so the same amount of information is contained in the spectrum (not in the power spectrum, that doesn't contain the phase, and hence, only half of the information)! Just convert your spectrum back to time domain.

However, if you just got the power spectrum plot of a signal, you'll find that many modulations have "characteristic" spectral shapes. With a bit of practice, you can tell binary GMSK from an 4-FSK pretty easily.

Now, in general, the problem of detecting a signal and its modulation is pretty hard – governments pay loads of money for well-working signal classificators, because in many cases, you're looking at a very weak $y(t)$ within a mixture of hundreds, and it's still hard to guess whether you're looking at an 16-QAM or an 8-PSK if you don't have timing information.

Thus: without a lot of additional info on what you're looking at, guessing the modulation class of $\alpha$ from $|Y(f)|$ is very hard. And it's harder to actually estimate $\alpha(t)$ from that.


Fourier transform of a damped cosine wave with a linear frequency chirp

For Fourier Transform of the LFM chirp portion - you use the Principle of Stationary Phase (POSP). The POSP essentially says the main contribution in the Fourier Integral comes from the portion of where the derivative of the phase is zero - it assumes that the integral of the oscillating components cancel themselves out.

Using the POSP - the Fourier transform of a LFM chirp is another LFM chirp in the frequency domain. The assumption is not very good in cases where you have a low time-bandwidth product.

For the complete signal, you have to options:

Evaluate the Fourier Transform of the magnitude portion and then convolve with the Fourier Transform of the LFM chirp pulse.
Keep the magnitude in while you are doing the POSP evaluation. You have to check on the assumptions the POSP make. If I recall correctly - it is usually assumes a slowly varying amplitude window.

I believe you can find the use of the POSP for LFM pulses in the signal analysis books by Papoulis, and also in "Digital Processing of Synthetic Aperture Radar Data" by Cummings and Wong.


Illumination Normalization Paper: Help understand equation

It is indeed a gradient but in particular a relative gradient or normalized gradient as the paper authors put it, they also give its definition:
$$
\nabla_rF(x,y) = \dfrac{\nabla F(x,y)}{|F(x,y)|+c}
$$, where $c$ is a constant to avoid division by zero (since $F(x,y)$ represents pixel intensities it is reasonable for a zero intensity value to occur, encoding a black pixel).
Because an image is represented in a discretized 2D spatial domain the gradient can be computed through finite differences around a neighborhood for each pixel's coordinates. A pixel in a rectilinear grid has at most 8 and at least 3 (for corner pixels) immediate neighbors, thus computing its relative differential involves taking 8 discrete differences, not including the difference with the center pixel itself.
Now the authors claim that this relative gradient is illumination insensitive because the image, in the context of lambertian illumination model, can be analyzed into two components: a reflectance and an illuminance component, with the latter varying slowly in a small neighborhood around a pixel and the first depending heavily on the characteristics of the face's surfaces.
$$
F(x,y) = R(x,y) \times I(x,y)
$$, where $R$ is the reflectance component and $I$ the illuminance component.
Combining the above, this leaves us with:
$$
WLD(F'(x,y)) = \arctan\left(\alpha\cdot\sum_{i \in A}\sum_{j \in A}\dfrac{F'(x,y)-F'(x-i\Delta x, y-j\Delta y)}{F'(x,y)+c}\right) \\
= \arctan\left(\alpha\cdot\sum_{i \in A}\sum_{j \in A}\dfrac{R(x,y)\times I(x,y)-R(x-i\Delta x, y-j\Delta y)\times I(x-i\Delta x, y-j\Delta y)}{R(x,y)\times I(x,y)+c}\right)
$$
Slow illuminance variations mean that we can locally approximate $I(x-i\Delta x, y-j\Delta y) \approx I(x,y) $.
Thus, by factoring out $I(x,y)$, ignoring $c$ and then adding a new constant $c'$:
$$
\arctan\left(\alpha\cdot\sum_{i \in A}\sum_{j \in A}\dfrac{R(x,y)\times I(x,y)-R(x-i\Delta x, y-j\Delta y)\times I(x-i\Delta x, y-j\Delta y)}{R(x,y)\times I(x,y)}\right)\\
\arctan\left(\alpha\cdot\sum_{i \in A}\sum_{j \in A}\dfrac{R(x,y)\times I(x,y)-R(x-i\Delta x, y-j\Delta y)\times I(x, y)}{R(x,y)\times I(x,y)}\right)\\
\arctan\left(\alpha\cdot\sum_{i \in A}\sum_{j \in A}\dfrac{R(x,y)-R(x-i\Delta x, y-j\Delta y)}{R(x,y)+c'}\right)
$$
which approximates the facial pixel intensity variations while being relatively invariant to intensity variations regarding the illuminance component of the image.
The final step of taking the arctangent of the product of the double sum and the scaling parameter $\alpha$, is used as a saturation step to avoid large values of the descriptors and to measure the direction of the variations.


FFT of SIN waves with different phase delays

You aren't looking at the frequency content of your sum of sinewave.  You are looking at the frequency content of a rectangular window on your sinewaves, and the window (FFT length) is shorter than the least common multiple of all the sinewave periods.

Signals that are not orthogonal within a window can cancel each other out partially or completely within that window.  What has happened is that you've chosen a window of a length such that each sinewave can get nearly completely cancelled out by the next higher and lower frequency pair of sinusoids of the "right" phase.  The ones on the ends are the only ones that aren't sandwiched, and thus aren't cancelled.


A very simple measure would be to compare each row in the image with the row above it, allowing some horizontal shift.

I've hacked together this simple algorithm in Mathematica:

Mean[MapThread[
  Function[{line1, line2},
   Min[Table[Norm[line1 - RotateLeft[line2, shift]], {shift, -5, 5}]]
   ], {s[[2 ;;]], s[[;; -2]]}]]

It simply takes each pair of adjacent rows, rotates one of the rows by -5..5 pixels, and takes the smallest euclidean distance. This yields one euclidean distance for each row pair. I simply take the mean (but depending on your actual data, a truncated mean or median might be more robust).

These are the results I get for artificially generated samples (Formula: Normalize(random noise * (1-factor) + signal * factor))

If I plot the result against the signal strength, the algorithm seems to measure the "meandering river signal strength" quite well:

EDIT: I forgot to normalize the input samples. Fixed that an uploaded new result images


Cochleagram vs STFT for CNN-based speech segregation

I was hoping someone would be able to give me an explanation
  hilighting some of the differences between cochleagrams and STFTs, and
  what their individual use-cases would be

STFT (short time fourier transform) is a technique that is used to obtain a time-frequency representation of audio signals. Spectrograms and phase vocoders are arguably the most important applications of the STFT.

The cochleagram, (which is really just a variant of the spectrogram), is also used to obtain a time-frequency representation of audio signals but has one major difference compared to STFTs (and spectrograms).

Regular spectrograms have a constant channel (frequency bin) bandwidth, which is given simply by $f_s / N$ (the sampling rate divided by the number of samples ). In contrast to that, cochleagrams' channels have different bandwidths, which are chosen to mimic human hearing  (the cochlea).
There're more than just one (empirical) equation to compute the cochleagram bands and one of them can be found here  ( https://en.wikipedia.org/wiki/Critical_band )

$$ERB(f) = 24.7 * (4.37 f / 1000 + 1) (eq 1)$$

In a cochleagram, the bandwidths grow increasingly larger with increasing frequency. Note also, that it's possible to create a cochleagram by using a normal STFT and then averaging the energy found in several frequency bins so that they match the bandwidth(s) given by eq #1. Alternatively, you can use gammatone filters, which are basically band-pass filters of different bandwidth.


The relationship [G/T] = [G]-[T] dB/K is incorrect.

G/T when in units are already given in dB/K is simply the gain in dB divided by the temperature in K.


How to find out the noise level of a signal

I just don't think you're thinking this through all the way, and using the toolbox in MATLAB before you understand the basics. SNR is signal-to-noise ratio and is defined as $\text{SNR}=\frac{\text{Signal Power}}{\text{Noise Power}}$. When you call snr(signal, noise) in MATLAB, all it is doing is calculating the signal power, mean(abs(signal).^2), and the noise power, mean(abs(noise).^2), taking the ratio of them and returning the result in decibels (dB).

In the comments you are asking about how to calculate the SNR of x (the noiseless signal)...well its dividing by zero because the noise power in a noiseless signal is zero. You'd be doing $\frac{\text{Signal Power}}{0}=\infty$, you can go ahead and try using the snr function by doing snr(signal, zeros(size(signal))) but I'm willing to bet it either returns Inf or gives an error.

Note: Read the MATLAB docs on snr(signal) and you'll realize it is probably not returning what you expect. Docs are here: https://uk.mathworks.com/help/signal/ref/snr.html


Estimating optimum sampling point for 16-QAM in MATLAB

An eye diagram can be created by just framing up your QAM signal (probably after matched filtering), where the frame length is the number of samples per QAM symbol. In MATLAB, this is very simple:

x = output_of_matched_filter;
N = number_of_samples_per_symbol;
eye = reshape(x,N,[]);
% "eye" now is a two-dimensional matrix, with each column containing one trace
% to plot on the eye diagram. you can then plot all of the traces with:
plot(eye);

One note: the call to reshape requires that the length of x is a multiple of N. If it isn't already, you'll need to slice the array to make it true:

nsyms = floor(length(x)/N);
eye = reshape(x(1:nsyms*N),N,nsyms);


I think if you added more details about your signal and the specifications of your hardware, that would be helpful to help you in that issue.
What you basically need is to multiply it by carrier signal and take the real part of the resulted signal, after that you need to do the inverse at the receiver end.  You might need also to upsample and use such filter equivalent to your system specifications.
Again, I think if you added more details about your system, we can also give you more information which may guide you.
Good luck


Should I use CV_HAAR_SCALE_IMAGE while using LBP CascadeClassifier?

The OpenCV documentation suggests using CASCADE_SCALE_IMAGE when performing a detectMultiScale with LBP features.

lollipops_cascade.detectMultiScale(frame_gray, lollipops,1.1, 2, 0|CASCADE_SCALE_IMAGE, Size(30, 30));

Even though, it works perfectly with  CV_HAAR_SCALE_IMAGE flag instead of CASCADE_SCALE_IMAGE.


If there is a function $H(e^{j\omega})$. Then will its conjugate be $H(e^{-j\omega})$?Why?

What you say only holds true for real signals. Let's prove it by the definition of Fourier transform:
$$
H(e^{j\omega}) = \mathcal{F}[h(t)] = \int_{-\infty}^{\infty} h(t) e^{-j\omega t} dt
$$
and $H^*(e^{j\omega})$ should be
$$
H^*(e^{j\omega}) = \int_{-\infty}^{\infty} h^*(t) e^{j\omega t} dt
$$
thus
$$
H^*(e^{-j\omega}) = \int_{-\infty}^{\infty} h^*(t) e^{-j\omega t} dt = \mathcal{F}[h^*(t)]
$$
If $h(n)$ is real, $\mathcal{F}[h^*(t)] = \mathcal{F}[h(t)] $, we can derive that
$$
H^*(e^{-j\omega}) = H(e^{j\omega})
$$
Take conjungate on both sides of the equation above, we have
$$
H(e^{-j\omega}) = H^*(e^{j\omega})
$$


The usage of the term "energy" in image processing context has historical reasons. It can be mapped to "brightness" or "intensity".

Imagine you have a (for matters of simplicity) greyscale image, that you want to transmit. Historically, that was done via some kind of analog modulation, where, simply speaking, dark "pixels", low intensity meant "low voltage" and bright "pixels", high intensity meant "high voltage". When looking at a CRT television, you could take the word "energy" quite literally, as an all black image will produce significantly lower energy consumption than an all white image.

The term "engergy minimisation" is usually employed, when trying to optimize some kind of algorithm or codec by comparing the difference of original signal and encoded-decoded signal. The codec would be perfect, if this difference was all zeros, so the target is "to minimise the energy" of the difference signal. Usage of "energy maximisation" is analog.

So, think of "energy" as "intensity" or "brightness" in the context of image processing.


Artifacts signal conversion from float to integer

The quantization error is dependent on the input signal. In your case the input signal is a sine wave and the quantization error is a combination of its 8 odd harmonics that fit into the available bandwidth. You can reduce the dependency of the quantization error to the input signal by adding triangular dither before quantizing: Add to each sample a uniformly distributed random number between 0 and the quantization step, and subtract from the result another uniformly distributed random number between 0 and the quantization step. Then quantize. This way you'll get a flat spectral noise floor (no distortion of the signal) and no modulation of the quantization noise by the signal.

For further reading on dither: Robert Alexander Wannamaker's PhD thesis "The Theory of Dithered Quantization", 2003.

Modifying your script to include dithering fixes the problem:

import numpy as np
import random
import matplotlib.pyplot as plt
from scipy.signal import get_window
from scipy.fftpack import fft

l = 128
f = 4
x = np.arange(l) / float(l)
max_int = np.iinfo(np.int16).max
sine = np.sin(2 * np.pi * f * x)
dither = [random.triangular(-1, 1, 0) for _ in xrange(l)]
sine_96 = (sine / sine.max() * (max_int - 1) + dither).astype(np.int16)

w = get_window('hann', l)
S = np.abs(fft(sine * w))
S_96 = np.abs(fft(sine_96 * w))

fig, ax = plt.subplots(2, 2, figsize=(8, 7))
ax[0, 0].plot(sine * w)
ax[0, 1].plot(sine_96 * w)

ax[1, 0].plot(20 * np.log10(S[:l // 2]))
ax[1, 1].plot(20 * np.log10(S_96[:l // 2]))
fig.show()

Figure 1. Top left: windowed input signal, top right: dithered, quantized, and windowed input signal, bottom left: frequency spectrum of the windowed input signal (vertical axis in dB scale), bottom right: frequency spectrum of the dithered, quantized, and windowed input signal (vertical axis in dB scale).

Let's plot the quantization residual (error) of the undithered case before windowing, and the frequency spectrum of the residual after windowing it:

sine_96 = (sine / sine.max() * max_int).astype(np.int16)
residual = sine_96 - (sine / sine.max() * max_int)
fig, ax = plt.subplots(1, 2, figsize=(8, 3))
ax[0].plot(residual)
R = np.abs(fft(residual * w))
ax[1].plot(20 * np.log10(R[:l // 2]))
fig.show()

Figure 2. Quantization residual and its spectrum (window function applied before transforming).

Quantization is equivalent to a memoryless waveshaper shaped like a staircase. For a periodic input signal a memoryless waveshaper can only produce harmonics of the input signal, because the resulting output signal will be identically periodic:

$$f(x + s) = f(x) \text{ for all } x \in \mathbf{R}\\
\Rightarrow g\big(f(x + s)\big) = g\big(f(x)\big) \text{ for all } x \in \mathbf{R}$$

Here the input signal is a sine wave with a period of 128/4 = 32 samples. If the period would not be an integer number of samples, then also other frequencies than the harmonics would be generated, because quantization operates in discrete time and non-integer periodicities don't exist in discrete time.


Fourier Transforms, Convolution, Cross-correlation: what is their physical unit exactly?

I agree that the units for Fourier transform are V/Hz. Most of my experience with the convolution integral dealt with a convolution of the impulse response of a linear network with the voltage input to the network. The impulse response for, say, a simple RC network will have the product RC in the denominator. Thus, the unit of impulse response is per second. So, the units of a convolution would be volts-seconds * per second = volts. For correlation, either auto or cross-, in the case of power signals (as opposed to energy signals), you should divide the integral by the period, T. If it is not periodic, still divide by T, but take the limit as T goes to infinity. Now the period, T, has units of seconds and you are integrating over time, which has units of seconds. The time in the numerator and denominator cancel, and you are left with volts squared, which gives power when divided by resistance (assuming that the volts are rms).


Relationship of Fourier Transform and Spatial Directivity Pattern

In your second link, the time and space dependence of a simple plane wave is given by

Ψ ( x , t ) = A sin ( ω t ± k x )

Just as a Fourier transform of the equation above can be taken from the time domain $t$ to the frequency domain $\omega$, so also can a Fourier transform be taken from the spatial $x$ domain to the wave vector $k$ domain. For the very simple one dimensional plane wave the wave vector reduces to the wavenumber  $k= \frac{2\pi}{\lambda}$ in that first equation.

In general both the position and wave vector are vector quantities and the dependence of a wave on both is given by

$$
e^{j\vec k\cdot\vec r}
$$

The aperture of a microphone array depends on the position of the microphones (a positional quantity), and so the Fourier transform depends on wave vector.

In Section 1.3.2 of your first link, they are doing exactly that - taking the Fourier transform of the aperture from position to wave vector space, which yields the pattern of the beam in terms of the wave vector.

In the example given in the other answer by @Peter K, think of a plane wave propagating away from the microphone array at an angle $\theta$ .The term inside the exponential, is the dot product of the wave vector and the spacing of the microphone array.


As mentioned in a comment, you can't use firls directly because this function only designs real-valued filters with (conjugate) symmetrical frequency responses. Your specification can only be met by a complex-valued filter because the frequency domain specification is not symmetrical.

In principle, you could use firls to design the real part and the imaginary part of the impulse response separately (by using the even and odd parts, respectively, of the frequency domain specification as the desired responses), but in that case you would get discontinuous desired responses, which would result in large approximation errors.

A better and simpler approach is to set up an overdetermined system of complex linear equations and use Matlab/Octave to solve them in a least squares sense. This is very straightforward if you use matrix commands:

N = 51;     % desired filter length
% frequency grid, desired frequency response, weighting function
f = [linspace(-1,-.18,164),linspace(-.1,.3,80),linspace(.38,1,124)];
d = [zeros(1,164),ones(1,80),zeros(1,124)].*exp(-1i*pi*f*(N-1)/2);
w = [2*ones(1,164),ones(1,80),2*ones(1,124)];
f = f(:); d = d(:); w = w(:);

% set up and solve overdetermined linear system
A = w(:,ones(1,N)) .* exp(-1i*pi*f*(0:N-1) );
h = A \ (w.*d);

The resulting filter is complex-valued and has a linear phase response, i.e., the real part of the impulse response is even, and the imaginary part is odd (see figure):

You can further decrease the approximation error by chosing a larger filter length.

I also have a corresponding Matlab/Octave function on GitHub: cfirls.m.


One option from the realms of adaptive filtering is the Least Mean Squares (LMS) filter depicted below:

The idea is you take the output of the unknown system, compare it with the output of your adaptive filter and minimize the difference by tweaking the filter coefficients, using a LMS algorithm. When the error $e(n)$ is zero (or more often, lower than a threshold you define) you found your unknown system's response.

There are countless resources where you can learn about LMS. In this case, given

the input vector: $$X(i) = [x(i) \space x(i-1) \ldots x(i-N+1)]^T$$ and the coefficients vector: $$W = [w_0 \space w_1 \ldots w_{N-1}]^T$$

The output of the adaptive filter is: $$\hat{y}(i)=X(i)^TW$$

Now define the cost function to be: $$e(i)^2 = \big(d(n)-\hat{y}(i)\big)^2$$
It can be shown that the gradient of the cost function with respect to the coefficients is:$$\hat{\nabla}(i) = \frac{\partial e(i)^2}{\partial W} = -2e(i)X(i)$$
Since the gradient points towards the maximum of the cost function you update the coefficients in the opposite direction: $$\hat{W}(i+1) = \hat{W}(i)-\mu\hat{\nabla}(i) = \hat{W}(i)+2\mu e(i)X(i)$$
where $\mu$ controls the adaptation speed. I've ignored a lot of details but this is the idea.

Note that the filter will produce a frequency response that is adapted to you input signal $x(n)$. It's not guaranteed for any other frequency. For example, in the response depicted below, my $x(n)$ is composed of 3 tones. You can see only at those particular frequencies does the filter (N=3,6,10) match the true response (H1(z)). So consider feeding it a wideband signal like white noise.

If you want to learn more you can check for example:

Wikipedia article

B. Widrow, Adaptive Signal Processing, Prentice-Hall, 1985


How do window size, sample rate influence FFT pitch estimation?

Use a Gaussian window - the Fourier transform of a Gaussian is a Gaussian
Log-scale the spectrum to emphasize peaks and turn the Gaussian peaks into parabolic peaks
Use parabolic interpolation to find the true peaks.

  Note that, as mentioned in §D.1, the Gaussian window transform magnitude is precisely a parabola on a dB scale. As a result, quadratic spectral peak interpolation is exact under the Gaussian window. Of course, we must somehow remove the infinitely long tails of the Gaussian window in practice, but this does not cause much deviation from a parabola, as shown in Fig.3.30.

https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks.html

I estimate 1000.000004 Hz for a 1000 Hz waveform this way: https://gist.github.com/255291#file_parabolic.py

If you're having trouble, plot the spectrum and use your eyes to see why it's not working.


What is leakage in frequency domain adaptive filters?

In adaptive filtering, leakage is a stabilization method which may be useful if the covariance matrix is close to singular (i.e. at least one of the eigenvalues is very small), or if there are finite-precision effects in the implementation of the adaptive filter. Leakage changes the update formula such that not only the mean squared error but also the norm of the filter taps is minimized. This prevents unbounded growth of the filter coefficients in cases of numerical ill-conditioning.

For you this simply means that you initially use a leakage factor of 1 (i.e. no leakage) if the FDAF works properly with your input signals. If you encounter coefficient drift (large fluctuation about the optimum solution), you can start by slightly decreasing the leakage factor until the coefficient fluctuation becomes sufficiently small. Note that leakage achieves stabilization at the expense of performance degradation because due to the changed update formula it introduces some bias in the filter taps.


Trouble Implementing Huggins Binaural Pitch

I just glanced at a paper on this subject, and it appears that you should be introducing an inter-aural 0-to-2PI phase shift over a very narrow frequency range, rather than offsetting one signal by T. I can imagine a variety of ways to create the phase-modified noise signal, included filtering a noise signal with a high-Q 2nd-order allpass filter, or by directly creating two noise signals by summing randomly-phased sinusoids, and directly imposing the desired phase(Left)-phase(right) condition by setting the phases of the 2 arrays appropriately.


How to calculate the noise power for a non-stationary noise?

If your signal is just $n(t)$ distributed as $N(\mu, \sigma^2)$, then the expected signal power is then:
$$
p_{\rm ensemble} = E\left [ \left| n(t) \right|^2 \right] = \mu^2 + \sigma^2
$$

For the time-varying case, I'd suggest that the answer is:
$$
p_{\rm time\ varying} = \mu^2(t) + \sigma^2(t)
$$


If the ROC is outside a circle in the complex $z$-plane ($|z|>a$), then the corresponding system is causal. If it is inside a circle ($|z|<a$), the system is anti-causal. If the ROC is a ring ($a<|z|<b$), the impulse response of the corresponding system is two-sided, i.e., the system is non-causal. In all cases, if the ROC contains the unit circle, the system is stable.

Consider the causal sequence $x[n]=a^nu[n]$ (where $u[n]$ is the discrete-time unit step sequence). The bilateral $\mathcal{Z}$-transform is defined by

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

Since $x[n]$ is causal we get

$$X(z)=\sum_{n=0}^{\infty}a^nz^{-n}=\sum_{n=0}^{\infty}\left(\frac{a}{z}\right)^n=\frac{1}{1-\left(\frac{a}{z}\right)}=\frac{z}{z-a},\qquad \left|\frac{a}{z}\right|<1\tag{2}$$

where I've used the formula for the geometric series. Note that the series only converges for $|a|<|z|$, i.e., outside a circle in the complex plane.

In a completely analogous manner we get for an anti-causal sequence such as $x[n]=a^nu[-n]$

$$X(z)=\sum_{n=-\infty}^{0}a^nz^{-n}=\sum_{n=0}^{\infty}a^{-n}z^n=\sum_{n=0}^{\infty}\left(\frac{z}{a}\right)^n=\frac{1}{1-\left(\frac{z}{a}\right)}=\frac{a}{a-z},\qquad \left|\frac{z}{a}\right|<1\tag{3}$$

which converges for $|z|<|a|$, i.e., inside a circle in the complex plane.

Finally, a two-sided sequence can be split up into a causal and an anti-causal sequence, with an ROC that is the overlap of the two individual ROCs, which, if it is not empty, is a ring in the complex plane.