FFT method input argument have to be $2^n$ ?

It is a common misconception that all FFT algorithms require their inputs to have a length that is a power of 2. The radix-2 Cooley-Tukey algorithm, the first to be widely known, has this limitation, but there are mixed-radix versions of the algorithm that don't require a power-of-2-sized input. In addition to Cooley-Tukey, there are a number of other FFT algorithms that have differing requirements on their input size; there are even techniques like Rader's algorithm that work on inputs whose sizes are prime.

Good, general purpose FFT implementations like MATLAB's fft() or FFTW support a wide range of algorithms and select one that is most appropriate for the specific situation. For more details on the specifics, go to the MATLAB documentation itself:

> doc fft

  Algorithms

  The FFT functions (fft, fft2, fftn, ifft, ifft2, ifftn) are based on a
  library called FFTW [3],[4]. To compute an N-point DFT when N is
  composite (that is, when N = N1N2), the FFTW library decomposes the
  problem using the Cooley-Tukey algorithm [1], which first computes N1
  transforms of size N2, and then computes N2 transforms of size N1. The
  decomposition is applied recursively to both the N1- and N2-point DFTs
  until the problem can be solved using one of several machine-generated
  fixed-size "codelets." The codelets in turn use several algorithms in
  combination, including a variation of Cooley-Tukey [5], a prime factor
  algorithm [6], and a split-radix algorithm [2]. The particular
  factorization of N is chosen heuristically.

  When N is a prime number, the FFTW library first decomposes an N-point
  problem into three (N – 1)-point problems using Rader's algorithm [7].
  It then uses the Cooley-Tukey decomposition described above to compute
  the (N – 1)-point DFTs.

  For most N, real-input DFTs require roughly half the computation time
  of complex-input DFTs. However, when N has large prime factors, there
  is little or no speed difference.

  The execution time for fft depends on the length of the transform. It
  is fastest for powers of two. It is almost as fast for lengths that
  have only small prime factors. It is typically several times slower
  for lengths that are prime or which have large prime factors.

    Note   You might be able to increase the speed of fft using the
    utility function fftw, which controls the optimization of the
    algorithm used to compute an FFT of a particular size and dimension.

  References

  [1] Cooley, J. W. and J. W. Tukey, "An Algorithm for the Machine
  Computation of the Complex Fourier Series,"Mathematics of Computation,
  Vol. 19, April 1965, pp. 297-301.

  [2] Duhamel, P. and M. Vetterli, "Fast Fourier Transforms: A Tutorial
  Review and a State of the Art," Signal Processing, Vol. 19, April
  1990, pp. 259-299.

  [3] FFTW (http://www.fftw.org)

  [4] Frigo, M. and S. G. Johnson, "FFTW: An Adaptive Software
  Architecture for the FFT,"Proceedings of the International Conference
  on Acoustics, Speech, and Signal Processing, Vol. 3, 1998, pp.
  1381-1384.

  [5] Oppenheim, A. V. and R. W. Schafer, Discrete-Time Signal
  Processing, Prentice-Hall, 1989, p. 611.

  [6] Oppenheim, A. V. and R. W. Schafer, Discrete-Time Signal
  Processing, Prentice-Hall, 1989, p. 619.

  [7] Rader, C. M., "Discrete Fourier Transforms when the Number of Data
  Samples Is Prime," Proceedings of the IEEE, Vol. 56, June 1968, pp.
  1107-1108.


Converting an arbitrary FIR filter into a linear-phase FIR with the same magnitude response is generally impossible. As mentioned in msm's answer, this conversion must be done with an allpass filter, which must be an IIR filter, and only in those rare cases where you get appropriate pole-zero cancellations will the total filter be FIR. In general you'll end up with an IIR filter.

Having said that, there are possibilities to convert a given FIR filter into a linear-phase FIR filter with approximately the same magnitude response. By increasing the length of the linear-phase FIR filter, the approximation can be made arbitrarily good.

The method is similar to the one you implemented, but the crucial difference is zero-padding. You need to use an FFT length that is much greater than the length of the given FIR filter.

Below is an octave example explaining how it works. I chose a non-linear phase filter of length $30$, and approximated it by a length $39$ linear-phase FIR filter:

% come up with some non-linear phase FIR filter
[b,a] = butter(5,.3);
h = impz(b,a);
h = h(1:30);    % length 30
h = h(:)';

% do the FFT trick
N = 1024;       % FFT length (even!)
H = fft(h,N);
hlp = real(ifft(abs(H)));
hlp = [hlp(N/2+2:N),hlp(1:N/2)];

% windowing
Nw = 39;        % choose window length (odd!) = filter length
w = ones(1,Nw); % rect. window (choose another if you wish)
w = w(:)';
Nw2 = (Nw-1)/2;
hlp = hlp(N/2-Nw2:N/2+Nw2) .* w;

% plot
Hlp = fft(hlp,N);
n = (0:N/2)*2/N;
subplot(2,1,1), stem(h), hold on, stem(hlp,'r'), hold off
legend('NLP FIR', 'LP FIR')
subplot(2,1,2),
plot(n,20*log10(abs(H(1:N/2+1))),n,20*log10(abs(Hlp(1:N/2+1))),'r')
legend('NLP FIR', 'LP FIR') ,grid on

The plot below shows the impulse responses of the two FIR filters, and the magnitude responses. The magnitude approximation could be further improved by choosing a longer filter length for the linear-phase FIR filter.


Goertzel is defined on complex numbers $I + jQ$, not real numbers. So, directly apply Goertzel to your baseband.

Your 1/128 scaling is unnecessary, since that's only a factor and doesn't contribute to the decision.

Also, your -127.4 (which almost certainly is a misunderstanding of some sorts) just kills the bias of unsigned integers coming out the device. Since that all ends up in the (not evaluated) DC bin of your DFT, you don't have to do that, either, and don't gain anything by doing so, either (numerical accuracy of floating point operations aside, but with any reasonable observation length, i.e. any reasonable bit rate, your accuracy cannot be degraded - after all, you start with just 8 bit of precision).

So, take these numbers, and feed them into your Goertzel, just converted to the data type your Goertzel implementation needs.

Another note: don't let yourself be fooled, computers are fast. 1 MS/s is not very much, and you need at most a full symbol-length FFT, i.e. a $\frac{10^6}{2400}=416\frac23\approx 417$-point FFT (yes, there's FFTs for non-power of two) will also take nearly no time. I've just benchmarked this, on my really elderly AMD CPU:

I only got 60,000 transforms per second, corresponding to 60,000 baud – that's only 25 as fast as you would actually need, so that's OK, but a bit disappointing.

Now, note that FFTs work better, the smaller the prime factors of their length is. Under that point of view, 1 MS/s is a terrible choice – it's not a multiple of 2400 /s. So, instead, use a much smaller sampling rate that is a multiple of 2.4 kHz. What about 240 kHz?
That would give you 100-FFTs, and of these, I get 1.3 million per second (you'd only need 2400 per second to follow your bits!).

So, the question really is why you came to the conclusion you needed 1MS/s for your 2400 bd signal in the first place. Unless it's absolutely unreasonably wide-spaced in spectrum, that's a waste of sampling rate. Your Dongle supports a wide range of sampling rates, among these multiples of your nominal baud rate – use such to keep your effort in range, or resample.

So, generally, regarding the question how to preprocess the data:

Use something that gives you samples that you can work with. I do like GNU Radio (that's what I used to build above flow graph), and it comes with all the tools you need. Even a weak personal computer-style CPU (like the ARM in a Raspberry Pi) will have plenty enough horsepower so that you can do a full FFT. Don't reinvent the wheel, though, but use existing implementations.

Of course, for detection of one of two potential frequencies, a Goertzel is actually more appropriate then an FFT. But if you're clever, you realize that actually you can get away with two Goertzels, one of which covers every frequency above 0, one every frequency below 0. That's not really any more efficient than a 3-point DFT, of which you just ignore the DC bin. So, I'd argue:

take your samples; don't scale or shift them, because that's just cosmetics (DC offset doesn't matter to a "is there more energy left or right in the spectrum, neither does scaling).
Like anyone who strives for acceptable SNR, you'd filter to only consider the part of the spectrum where there might be your signal. Use a real-tapped (ie. spectrum-symmetrical) band pass filter to

suppress the DC component we didn't suppress, alongside with LO leakage
select your bands of interest from noise

Notice that above filter might optimally be (depends on whether you only consider an AWGN channel) the conjugate pulse shaping filter, shifted up with a cosine
Decimate the output of that filter. You can decimate very drastically here, just make sure if aliasing occurs, it always happens within the same band where the original energy came from. If you're not sure you can achieve that, only decimate to Nyquist rate.
Circularly shift in frequency (by multiplying with a $e^{j\frac{2\pi}{x}t} $) so that one of your tones now is at DC, and the other is not.
Your 2 Goertzels now collapse to the full 2-DFT, which really collapses to

First bin: Moving average of samples
Second bin: Difference of successive samples


The soundsc command will playback $f_s$ samples per second. Playing back at a faster rate than you recorded will increase the pitch of the speech. Playing back at a lower sample rate will reduce the pitch.

If you prefer to think in the frequency domain; all the information is scaled on the frequency axis with relation to the change in $f_s$. Reducing fs squashes the information on the frequency axis and increasing $f_s$ will stretch the information on the frequency axis.


Problem with scipy.signal.butter for designing filter bank

Although Im not sure how I can display the frequency response using this method, I have found that it is workable to use second-order section method for filtering. So say we have a signal x with sampling frequency fs_real that needs to filtered using a constant percentage bandwidth (in this case third octaves) then the corresponding power spectrum will be given by the get_cpb function:

def butter_bandpass_sos(lowcut, highcut, fs, order=9):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    sos = butter(order, [low, high], btype='band', output='sos')
    return sos

def get_cpb(x,fs_real):
    """
    Calculate third octave spectrum from pressure time series.
    An octave is the interval between two frequencies having
    a ratio of 2:1. Hence for 1/3 octave
    fh = fc*(sqrt(2)^(1/3))
    fl = fc/(sqrt(2)^(1/3))
    """

    cpb = []

    for fc in F0:
        fh = fc*(2**(1/6.0))
        fl = fc/(2**(1/6.0))
        sos = butter_bandpass_sos(fl,fh,fs_real)
        y = sosfilt(sos,x)
        cpb.append( 20*np.log10(np.std(y)/2e-5) )

    return F0, cpb


How to flatten the image of a label on a food jar?

A similar question was asked on Mathematica.Stackexchange. My answer over there evolved and got quite long in the end, so I'll summarize the algorithm here.

Abstract

The basic idea is:

Find the label.
Find the borders of the label
Find a mapping that maps image coordinates to cylinder coordinates so that it maps the pixels along the top border of the label to ([anything] / 0), the pixels along the right border to (1 / [anything]) and so on.
Transform the image using this mapping

The algorithm only works for images where:

the label is brighter than the background (this is needed for the label detection)
the label is rectangular (this is used to measure the quality of a mapping)
the jar is (almost) vertical (this is used to keep the mapping function simple)
the jar is cylindrical (this is used to keep the mapping function simple)

However, the algorithm is modular. At least in principle, you could write your own label detection that does not require a dark background, or you could write your own quality measurement function that can cope with elliptical or octagonal labels.

Results

These images were processed fully automatically, i.e. the algorithm takes the source image, works for a few seconds, then shows the mapping (left) and the un-distorted image (right):

The next images were processed with a modified version of the algorithm, were the user selects the left and right borders of the jar (not the label), because the curvature of the label cannot be estimated from the image in a frontal shot (i.e. the fully automatic algorithm would return images that are slightly distorted):

Implementation:

1. Find the label

The label is bright in front of a dark background, so I can find it easily using binarization:

src = Import["http://i.stack.imgur.com/rfNu7.png"];
binary = FillingTransform[DeleteBorderComponents[Binarize[src]]]

I simply pick the largest connected component and assume that's the label:

labelMask = Image[SortBy[ComponentMeasurements[binary, {"Area", "Mask"}][[All, 2]], First][[-1, 2]]]

2. Find the borders of the label

Next step: find the top/bottom/left/right borders using simple derivative convolution masks:

topBorder = DeleteSmallComponents[ImageConvolve[labelMask, {{1}, {-1}}]];
bottomBorder = DeleteSmallComponents[ImageConvolve[labelMask, {{-1}, {1}}]];
leftBorder = DeleteSmallComponents[ImageConvolve[labelMask, {{1, -1}}]];
rightBorder = DeleteSmallComponents[ImageConvolve[labelMask, {{-1, 1}}]];

This is a little helper function that finds all white pixels in one of these four images and converts the indices to coordinates (Position returns indices, and indices are 1-based {y,x}-tuples, where y=1 is at the top of the image. But all the image processing functions expect coordinates, which are 0-based {x,y}-tuples, where y=0 is the bottom of the image):

{w, h} = ImageDimensions[topBorder];
maskToPoints = Function[mask, {#[[2]]-1, h - #[[1]]+1} & /@ Position[ImageData[mask], 1.]];

3. Find a mapping from image to cylinder coordinates

Now I have four separate lists of coordinates of the top, bottom, left, right borders of the label. I define a mapping from image coordinates to cylinder coordinates:

arcSinSeries = Normal[Series[ArcSin[\[Alpha]], {\[Alpha], 0, 10}]]
Clear[mapping];
mapping[{x_, y_}] :=
   {
    c1 + c2*(arcSinSeries /. \[Alpha] -> (x - cx)/r) + c3*y + c4*x*y,
    top + y*height + tilt1*Sqrt[Clip[r^2 - (x - cx)^2, {0.01, \[Infinity]}]] + tilt2*y*Sqrt[Clip[r^2 - (x - cx)^2, {0.01, \[Infinity]}]]
   }

This is a cylindrical mapping, that maps X/Y-coordinates in the source image to cylindrical coordinates. The mapping has 10 degrees of freedom for height/radius/center/perspective/tilt. I used the Taylor series to approximate the arc sine, because I couldn't get the optimization working with ArcSin directly. The Clip calls are my ad-hoc attempt to prevent complex numbers during the optimization. There's a trade-off here: On the one hand, the function should be as close to an exact cylindrical mapping as possible, to give the lowest possible distortion. On the other hand, if it's to complicated, it gets much harder to find optimal values for the degrees of freedom automatically. (The nice thing about doing image processing with Mathematica is that you can play around with mathematical models like this very easily, introduce additional terms for different distortions and use the same optimization functions to get final results. I've never been able to do anything like that using OpenCV or Matlab. But I never tried the symbolic toolbox for Matlab, maybe that makes it more useful.)

Next I define an "error function" that measures the quality of a image -> cylinder coordinate mapping. It's just the sum of squared errors for the border pixels:

errorFunction =
  Flatten[{
    (mapping[#][[1]])^2 & /@ maskToPoints[leftBorder],
    (mapping[#][[1]] - 1)^2 & /@ maskToPoints[rightBorder],
    (mapping[#][[2]] - 1)^2 & /@ maskToPoints[topBorder],
    (mapping[#][[2]])^2 & /@ maskToPoints[bottomBorder]
    }];

This error function measures the "quality" of a mapping: It's lowest if the points on the left border are mapped to (0 / [anything]), pixels on the top border are mapped to ([anything] / 0) and so on.

Now I can tell Mathematica to find coefficients that minimize this error function. I can make "educated guesses" about some of the coefficients (e.g. the radius and center of the jar in the image). I use these as starting points of the optimization:

leftMean = Mean[maskToPoints[leftBorder]][[1]];
rightMean = Mean[maskToPoints[rightBorder]][[1]];
topMean = Mean[maskToPoints[topBorder]][[2]];
bottomMean = Mean[maskToPoints[bottomBorder]][[2]];
solution =
 FindMinimum[
   Total[errorFunction],
    {{c1, 0}, {c2, rightMean - leftMean}, {c3, 0}, {c4, 0},
     {cx, (leftMean + rightMean)/2},
     {top, topMean},
     {r, rightMean - leftMean},
     {height, bottomMean - topMean},
     {tilt1, 0}, {tilt2, 0}}][[2]]

FindMinimum finds values for the 10 degrees of freedom of my mapping function that minimize the error function. Combine the generic mapping and this solution and I get a mapping from X/Y image coordinates, that fits the label area. I can visualize this mapping using Mathematica's ContourPlot function:

Show[src,
 ContourPlot[mapping[{x, y}][[1]] /. solution, {x, 0, w}, {y, 0, h},
  ContourShading -> None, ContourStyle -> Red,
  Contours -> Range[0, 1, 0.1],
  RegionFunction -> Function[{x, y}, 0 <= (mapping[{x, y}][[2]] /. solution) <= 1]],
 ContourPlot[mapping[{x, y}][[2]] /. solution, {x, 0, w}, {y, 0, h},
  ContourShading -> None, ContourStyle -> Red,
  Contours -> Range[0, 1, 0.2],
  RegionFunction -> Function[{x, y}, 0 <= (mapping[{x, y}][[1]] /. solution) <= 1]]]

4. Transform the image

Finally, I use Mathematica's ImageForwardTransform function to distort the image according to this mapping:

ImageForwardTransformation[src, mapping[#] /. solution &, {400, 300}, DataRange -> Full, PlotRange -> {{0, 1}, {0, 1}}]

That gives the results as shown above.

Manually assisted version

The algorithm above is full-automatic. No adjustments required. It works reasonably well as long as the picture is taken from above or below. But if it's a frontal shot, the radius of the jar can not be estimated from the shape of the label. In these cases, I get much better results if I let the user enter the left/right borders of the jar manually, and set the corresponding degrees of freedom in the mapping explicitly.

This code lets the user select the left/right borders:

LocatorPane[Dynamic[{{xLeft, y1}, {xRight, y2}}],
 Dynamic[Show[src,
   Graphics[{Red, Line[{{xLeft, 0}, {xLeft, h}}],
     Line[{{xRight, 0}, {xRight, h}}]}]]]]

This is the alternative optimization code, where the center&radius are given explicitly.

manualAdjustments = {cx -> (xLeft + xRight)/2, r -> (xRight - xLeft)/2};
solution =
  FindMinimum[
   Total[minimize /. manualAdjustments],
    {{c1, 0}, {c2, rightMean - leftMean}, {c3, 0}, {c4, 0},
     {top, topMean},
     {height, bottomMean - topMean},
     {tilt1, 0}, {tilt2, 0}}][[2]]
solution = Join[solution, manualAdjustments]


How to Fourier Transform (basic understanding)

If you want to plot the magnitude of your data you need to calculate first the absolute value as the output of the fft is complex.

z = fft(y);
z_mag=abs(z); % you could also log-tranform the data for db-scale

To get the right x-Axis you need however the sampling rate. Lets assume that the sampling rate at which the data is recored is f_sampling(some number);

Fs=f_sampling;
n=length(z_mag)
freqHz = (0:1:n-1)*Fs/N;

Keep in mind that this is a double-sided spectrum.


Unable to obtain correct correct fourier transform for homework

If $x(t)$ is the (non-periodic) function in the interval $[-2,2]$, then the periodic function $x_p(t)$ is given by

$$x_p(t)=x(t)\star \sum_{n=-\infty}^{\infty}\delta(t-6n)\tag{1}$$

where $\star$ denotes convolution, and, consequently, its Fourier transform is

$$X_p(f)=X(f)\cdot\frac16\sum_{k=-\infty}^{\infty}\delta\left(f-\frac{k}{6}\right)=
\frac16\sum_{k=-\infty}^{\infty}X\left(\frac{k}{6}\right)\delta\left(f-\frac{k}{6}\right)\tag{2}$$

where $X(f)$ is the Fourier transform of $x(t)$. It remains to find $X(f)$. The function $x(t)$ can be written as

$$x(t)=\sin\left(\frac{\pi t}{2}\right)\cdot\left[\text{rect}\left(\frac{t-1}{2}\right)-\text{rect}\left(\frac{t+1}{2}\right)\right]\tag{3}$$

With the Fourier transform pairs

$$\begin{align}\sin\left(\frac{\pi t}{2}\right)&\Longleftrightarrow \frac{1}{2j}\left[\delta\left(f-\frac14\right)-\delta\left(f+\frac14\right)\right]\\\text{rect}\left(\frac{t\pm 1}{2}\right)&\Longleftrightarrow 2e^{\pm j2\pi f}\text{sinc}(2f)\end{align}\tag{4}$$

we obtain from $(3)$

$$\begin{align}X(f)&=\frac{1}{2j}\left[\delta\left(f-\frac14\right)-\delta\left(f+\frac14\right)\right]\star 2\;\text{sinc}(2f)(e^{-j2\pi f}-e^{j2\pi f})\\&=2\left[\delta\left(f+\frac14\right)-\delta\left(f-\frac14\right)\right]\star\text{sinc}(2f)\sin(2\pi f)\\
&=2\left[\text{sinc}\left(2f+\frac12\right)\sin\left(2\pi f+\frac{\pi}{2}\right)-\text{sinc}\left(2f-\frac12\right)\sin\left(2\pi f-\frac{\pi}{2}\right)\right]\end{align}\tag{5}$$

The final expression for $X_p(f)$ results from combining $(5)$ and $(2)$.


How to quantify amplitude of oscillations in marginally stable systems?

Is amplitude of oscillation a (only) function of initial condition of the system?

No.  For a system with a transfer function of the form $H(s) = \frac{\cdots}{s^2 + \omega_0^2}$, the amplitude of oscillation will increase any time it is excited by a signal that has a component at $x(t) = \cos \omega_0 t + \phi$.  Even if you're not exciting it intentionally, random noise will excite that pole.  Basically, the output will be of the form $a(t) \cos \omega_0 t + b(t) \sin \omega_0 t$, where $a(t)$ and $b(t)$ will be Wiener (random-walk) processes.  Such processes have variances that tend to infinity as $t \to \infty$.

How to intuitively quantify amplitude of the oscillations for marginally stable system?

"Big and growing, boss, and I don't think they'll get smaller!"
Possibly by how fast they grow.  But grow they will.

In some applications, oscillations with small amplitude might be acceptable.

If you have a system that's exhibiting persistently small oscillations, then you're seeing the signs of a nonlinear phenomenon called a "limit cycle".  Basically, there's some oscillator in there that's either inherently small-valued (like the oscillation you may see around the least significant bit of a DAC or ADC), or there's a big oscillation that's mostly not getting to the output.  Either way, if it persists at some amplitude, there's some nonlinear process that's keeping it that way.


Make use of perceptual hashes. It is very fast to compute and is very lightweight, both in terms of memory and cpu consumption. They are represented by simple long integers and can be indexed using many types of data structures such as VP Trees:

http://www.phash.org/

If that doesn't work, you can extract SURF features, quantize them into visual words using a vocabulary tree or BOW and index these with a database structure such as inverted file index. Check out OpenCV. It has many of the tools I mentioned here.


A DFT includes an implicit finite length window on any function (which has an infinite domain in the case of a sinusoid).  The effect of multiplication by that window in the time domain on the DFT result is a convolution with the transform of that window.  In the case of a rectangular window, that looks like something very close to a Sinc function (except near DC or Fs/2), or sin(w)/w (which has lots of "ups and downs".)


How to choose the most appropriate color space for a computer vision application?

Computer vision is a huge topic and your particular example (statistical color transfer) is too narrow to provide guidance on the remaning parts.

That being said, classically most algorithms have worked on black and white "color" space, which was either obtained directly from the sensor or was later converted to it. The reasons behind is: 1-It's much simpler, 2- It provides remarkably good results for most applications unless color is explicitly required.

Color is, however, inevitably required for certain machine vision applications. Then comes the problem of selecting the most suitable color space. RGB is the fundamental color space because those R,G,B signals are physically captured by the sensor. (This also came from human color perception and additive color formation based on it)

Correlation in RGB channels should be understood clearly. In general given an R value you cannot say anything about G or B. But given a brightness value and hue, then RGB will (may) be correlated.

YUV is the most useful color space used in TV broadcast technology, image and video coding applications (codecs, cameras, camcorders) and many technical applications where color is necessary. For some reasons YUV was derived based on special requirements of color tv broadcasting over existing black and white tvs with compatibility in mind. Y channel carries the brightness of the scene and typically used to replace the black and white image. U and V channels carry color value and saturation information.

HSV is a more complicated color space to derive from RGB. It can provide more accurate color discrimination than YUV. The HSV is definining image aspects like Hue, Saturation and Brightness which are fundamentally uncorrelated (or less correlated compared to RGB)

For technical reasons other color spaces also exist but they do not provide any advantages for signal processing and especially for machine vision.


Implementation of Wikipedia Equation for the DFT

You have a bug in ft2.  You are incrementing i, and freq together.  That's not how you want your summation to work.  I messed around with fixing it, but it got messy.  I decided to rewrite it from a discrete perspective instead of trying to use the continuous terminology.  In the DFT, the sampling rate is irrelevant.  What matters is how many samples are used (N).  The bin numbers (k) then correspond to frequency in units of cycles per frame.  I tried to keep you code as intact as possible so it would remain easily comprehensible to you.  I also unfurled the DFT calculation loops to hopefully reveal their nature a little bit better.

Hope this helps.

Ced

import numpy as np
import matplotlib.pyplot as plt

def ft(t, s, fs):
    freq_step = fs / len(s)
    freqs = np.arange(0, fs/2, freq_step)
    S = []
    for freq in freqs:
        real = np.sum(s * np.cos(2*np.pi*freq * t))
        compl = np.sum(- s * np.sin(2*np.pi*freq * t))
        tmpsum = (real**2 + compl**2) ** 0.5
        S.append(tmpsum)
    return S, freqs

def ft3(s, N):  # More efficient form of wikipedia equation

    S = []

    slice  = 0.0
    sliver = 2*np.pi/float(N)

    for k in range(N/2):

        sum_real = 0.0
        sum_imag = 0.0
        angle = 0.0
        for n in range(N):
            sum_real +=  s[n] * np.cos(angle)
            sum_imag += -s[n] * np.sin(angle)
            angle += slice

        slice += sliver
        tmpsum = (sum_real**2 + sum_imag**2) ** 0.5
        S.append(tmpsum)

    return S

def ft4(s, N):  # Using wikipedia equation

    S = []

    for k in range(N/2):

        sum_real = 0.0
        sum_imag = 0.0
        for n in range(N):
            sum_real +=  s[n] * np.cos(2*np.pi*k*n/float(N))
            sum_imag += -s[n] * np.sin(2*np.pi*k*n/float(N))

        tmpsum = (sum_real**2 + sum_imag**2) ** 0.5
        S.append(tmpsum)

    return S

def ft5(s, N):  # Roots of Unity Weighted Sum

    sliver = 2 * np.pi / float( N )

    root_real = np.zeros( N )
    root_imag = np.zeros( N )

    angle = 0.0
    for r in range(N):
        root_real[r] =  np.cos( angle )
        root_imag[r] = -np.sin( angle )
        angle += sliver

    S = []

    for k in range( N/2 ):

        sum_real = 0.0
        sum_imag = 0.0
        r = 0

        for n in range( N ):
            sum_real += s[n] * root_real[r]
            sum_imag += s[n] * root_imag[r]
            r += k
            if r >= N : r -= N

        tmpsum = np.sqrt( sum_real*sum_real + sum_imag*sum_imag )
        S.append( tmpsum )

    return S

def main():

    N  = 200
    fs = 100.0

    time_step = 1.0 / fs
    t = np.arange(0, N * time_step, time_step)

    f = 5.0
    y = np.sin(2*np.pi*f*t) + np.cos(2*np.pi*f*2*t)

    fig = plt.figure()
    ax = fig.add_subplot(311)
    ax.set_title('Signal in time domain')
    ax.set_xlabel('t')
    ax.plot(t, y)

    S, freqs = ft(t, y, fs)

    ax = fig.add_subplot(312)
    ax.set_xticks(np.arange(0, freqs[-1], 2))
    ax.set_title('Time using equation')
    ax.set_xlabel('frequency')
    ax.plot(freqs, S)

    S = ft3(y, N)
    ax = fig.add_subplot(313)
    ax.set_title('Using Wiki equation')
    ax.set_xlabel('frequency')
    ax.set_xticks(np.arange(0, freqs[-1], 2))
    print len(S), len(freqs)
    ax.plot(freqs, S)

    plt.tight_layout()
    plt.show()

main()


What is the Difference between Texture and Image?

An image (2D) is the projection of a natural scene or man made drawing on to a plane. The output of an optical sensor when displayed is an image.

Texture is a picture created by repeating a basic element called 'texel'. The texel can be a natural image or an artificial one.

The edge detection algorithms works on matrices. So if it is a normal image or texture image, you can do in the same way.

In texture images, recognizing the texel makes the rest of processing easy. Because you need to work on a smaller sized matrix.


Power spectral density: Why are these two methods equal?

Regarding what we discussed in the comments, here is an explanation of the different concepts I was referring to:

The concrete usage of terms may vary, but my interpretation is this: For a stochastic process $x(t)$, the autocorrelation function is defined as $$r(\tau) = \mathbb{E}\{x(t) \cdot x(t+\tau)\}.$$ Note that the process needs to be stationary, otherwise the autocorrelation function depends on $t$ as well. Its power spectral density is given by $$R(f) = \mathbb{E}\{|X(f)|^2\}.$$ These two quantities are linked to each other via the Wiener-Khinchin-Theorem which shows that $$R(f) = \int_{-\infty}^{\infty} r(\tau) {\rm e}^{-\jmath 2\pi f \tau} {\rm d} \tau,$$ i.e., they are Fourier transform pairs.

Since the expectation is hard to implement in practice (we typically do not have an infinite ensemble available), we need to resort to other methods. When the process is ergodic, ensemble average can be replaced by temporal average and we may, e.g., estimate the autocorrelation function from a window $T$ via $$r_T(\tau) = \frac{1}{T} \int_T x(t) \cdot x(t+\tau){\rm d}t.$$ Note that $r(\tau) = \lim_{T\rightarrow \infty} r_T(\tau)$. On the other hand, considering using $R_T(f) = |X_T(f)|^2$ as an estimate for the power spectral density (referred to as the periodigram), where $X_T(f)$ is the spectum of a $T$-windowed signal leads to an estimate that does not converge since $S(f) \neq \lim_{T\rightarrow \infty} S_T(f)$.

Now, I'm not sure what exactly is the question. Is it the explanation of the Wiener-Khinchin-Theorem for the stochastic signals? Or is it the deterministic equivalent as, e.g., discussed in this answer?

*edit: Regarding your clarification that you posted as a reply to your question(?): My feeling is this question has little to do with ACF and PSD being Fourier transform pairs, instead you seem to be worried how to normalize the result of a Discrete Fourier transform (carried out via an FFT) such that its units make sense physically. It's probably worth asking this in a separate question, but the short version is here. There are two things to be considered:

There are different conventions how to normalize DFT and IDFT. Let $\mathbf{F}$ be a DFT matrix with elements ${\rm e}^{-\jmath \frac{mn}{N}}$, then you can either (i) define DFT via $\mathbf{D} = \mathbf{F} \cdot \mathbf{d}$ in which case you need to define the IDFT via $\mathbf{d} = {\color{red}{ \frac{1}{N} }} \mathbf{F}^H \cdot \mathbf{D}$ or (ii) define the DFT via $\mathbf{D} = {\color{red}{ \frac{1}{N} }}\mathbf{F} \cdot \mathbf{d}$ which leads to the IDFT $\mathbf{d} =  \mathbf{F}^H \cdot \mathbf{D}$ or alternatively (iii) you let the DFT be $\mathbf{D} = {\color{red}{ \frac{1}{\sqrt{N}} }}\mathbf{F} \cdot \mathbf{d}$ which leads to the IDFT $\mathbf{D} = {\color{red}{ \frac{1}{\sqrt{N}} }} \mathbf{F}^H \cdot \mathbf{d}$. Different conventions are in place, Matlab uses the first by default, whereas Maple will prefer the third.
The Fourier transform changes the units: the ACF has a unit of power (say, Watt), the Fourier transform is an integral over time, hence it adds seconds, so that the unit is Watt seconds or better Watt/Hz. It's a density anyways (how do Watts distribute over Hertz?). If you want to approximate a Fourier transform via a Discrete Fourier transform, you have to be super careful, since the DFT/IDFT does not change the units. If your input to the DFT is a sequence with unit, say, Watt, its DFT coefficients will have the same unit. If you want to mirror what happens in a Fourier transform, you need to multiply the result with the sampling period of the samples in time domain. In other words: If you have a function function $u(t)$ with spectrum $U(f)$ and you want to approximate $U(f)$ with a DFT, you need to consider a finite window (say, length $T$) and a finite sampling interval (say, $t_0 = T/N$). Then, ${\color{red}{t_0}}\cdot D[\mu]$ is an approximation for $U(\mu f_0)$ where $f_0 = 1/T$. It's an approximation since sampling as well as truncation can cause errors (one of them will for sure).

I think the factor you were looking for can very likely be a combination of these two effects, your $1/(f_s N)$ is the same as my $t_0/N$.


Artificial Neural Network Preprocessing in Real-Time Applications

The answer i give is probably partly help full.
What you can do is when you have a pre stored STD and MEAN(N samples long). While running the real time proces update the stored MEAN.
Create an array of you mean value. During the loop add every cycle a sample to your array and remove the oldest one. You create a kind of circular pattern. After N samples your MEAN is only the MEAN from your data.
Do the same with the STD.

Maybe you can create a different MEAN and STD from different types of data so the error is less worse

Good Luck


Presented are two vectors of length $l$ with a displacement of $k>0$ between them

$A_i=a_i+ω_i$

$B_i=b_i+ω_i$

where $ω_i\sim N(0,σ_n^2)$ and $a_i\sim N(0,σ_s^2)$. The displacement means that for $k<i<l-k$, $b_{i+k}=a_i$

Our estimation of the circular correlation, without zero-padding, takes the form of

$C_k=\frac{∑_{i=1}^NA_i B_{mod(i-k+m,m)}}{\sqrt{∑_{i=1}^NA_i^2∑_{i=1}^NB_i^2}} \frac{l}{l-k}$

Following derivation

$C_k=\frac{∑_{i=1}^N (a_i+ω_i)(b_{mod(i-k+m,m)}+ω_{mod(i-k+m,m)})}{\sqrt{∑_{i=1}^N(a_i+ω_i)^2∑_{i=1}^N(b_i+ω_i)^2}} \frac{l}{l-k}$

Expected value

The expected value for square of one sample is

$E[(a_i+ω_i )^2 ]=E[a_i^2+a_i ω_i+ω_i^2 ]=σ_s^2+σ_n^2$

Expected value for the circular correlation at the true peak location
$E[C_k]=\frac{σ_s^2}{σ_s^2+σ_n^2}$

which means that we use biased estimation of the true peak correlation.

Correlation value not at the true peak location is
$E[C_c]=0$

Variance

General expression for the variance of multiplication of two samples is given by
$var((a_i+ω_i)(b_i+ω_j))=var(a_i b_i+ω_i b_i+a_i ω_i+ω_i ω_j )$

For the correlated samples

$var(a_ib_i+ω_ib_i+a_iω_i+ω_iω_j)=var(2a_iω_i+ω_iω_j)=4σ_s^2 σ_n^2+σ_n^4$

For the uncorrelated samples

$var(a_ib_i+ω_ib_i+a_iω_i+ω_iω_j)=var(a_ib_i+2a_iω_i+ω_iω_j)=σ_s^4+4σ_s^2σ_n^2+σ_n^4$

Lets neglect the variance of the dynamometer. This leads to the variance of the true peak correlation value

$var(C_k )≈\left(k(σ_s^4+4σ_s^2 σ_n^2+σ_n^4 )+(l-k)(4σ_s^2 σ_n^2+σ_n^4)\right)\left(\frac{l}{l-k}\frac{1}{l(σ_s^2+σ_n^2)}\right)^2$

simplification leads to

$var(C_k)=\frac{kσ_s^4+l(4σ_s^2 σ_n^2+σ_n^4)}{(l-k)^2(σ_s^2+σ_n^2 )^2}$

and for a location which is not the true peak v≠k

$var(C_v)≈l(σ_s^4+4σ_s^2σ_n^2+σ_n^4)\left(\frac{l}{l-v}\frac{1}{l(σ_s^2+σ_n^2)}\right)^2$

here simplification leads to

$var(C_v)=\frac{l(σ_s^4+4σ_s^2σ_n^2+σ_n^4)}{(l-v)^2(σ_s^2+σ_n^2)^2}$

Therefore, the true peak value

$C_k\sim\left(\frac{σ_s^2}{σ_s^2+σ_n^2},\frac{kσ_s^4+l(4σ_s^2 σ_n^2+σ_n^4)}{(l-k)^2(σ_s^2+σ_n^2 )^2}\right)$

and the false peak value

$C_v\sim\left(0,\frac{l(σ_s^4+4σ_s^2σ_n^2+σ_n^4)}{(l-v)^2(σ_s^2+σ_n^2)^2}\right)$


How does a FIR equiripple filter behave close to DC?

If you want to know the exact DC gain of your filter take the sum of the coefficients. This is true for an equiripple filter as well as any other FIR filter.  There is no such thing as ripple at DC, because "ripple" is a term used to describe a magnitude variation over some range of frequencies.

If what you are looking for is a flat passband then equiripple can fit the bill.  When using the fdatool in Matlab, notice that the equiripple FIR filter allows you to specify a parameter Apass.  This paramter lets you adjust the maximum allowable amplitude of the ripples in the passband as shown below.  By default it is set to 1, which means that you will have a maximum ripple of 1 dB in the passband.

If you are not constrained to FIR filters, then you might want to consider an IIR filter, depending on your needs.  You can usually achieve the desired result with a lower filter order when using IIR.  This means less computation and less delay.


Confusion in basics of Laplace Transform

The unilateral Laplace transform is used for analyzing causal linear time-invariant systems, which have an impulse response $h(t)$ that is zero for $t<0$. The unilateral Laplace transform can be used to solve initial value problems, due to the correspondence

$$x'(t)\Longleftrightarrow sX(s)-x(0)$$

where $x(0)$ is a given initial value for the function $x(t)$. Note that for the bilateral Laplace transform the equivalent correspondence is simply $x'(t)\Longleftrightarrow sX(s)$.

Concerning the signal $x(t)=e^{-at}u(t)$, note that its Laplace transform is

$$X(s)=\int_0^{\infty}e^{-(a+s)t}dt$$

This integral only converges if the exponential decays, which results in the condition $\text{Re}(a+s)>0$. This condition defines the region of convergence (ROC) of the Laplace transform. $X(s)$ only exists for values of the complex variable $s$ satisfying $\text{Re}(s)>-\text{Re}(a)$.


Your image is M rows x N cols, and a 2D FFT requires 1D FFTs on all rows followed by 1D column-wise FFTs on the result.

So you have:

$M N \log N + N M \log M = MN (\log M + \log N) = MN \log MN$

operations in total.

The DFT case can be derived using a similar method - this is left as an exercise for the reader.


Find the impulse response of an LTI system?

For the case a-)
I assume you are trained enough in DSP to see that :
$$ \mathcal{ F } \{ \delta[n-d] \} = e^{-j \omega d } $$ for all integer $d$.

Hence given a frequency response of $H(\omega) = e^{-j3\omega }$ it's apparent that the corresponding impulse response is $$h_a[n] = \delta[n-3] $$

The problem in case b-) is in the fact that it suggests a non-integer amount of shift $\alpha$ of the unit impulse $\delta[n]$ so as to obtain $\delta[n-\alpha]$ which would then have a corresponding frequency response $H_b(\omega) = e^{ -j \omega \alpha }$. But this makes no sense in the domain of discrete-time sequences which cannot be shifted by non-integer amounts which might force you
to argue that the corresponding impulse response does not even exist.

The solution requires an investigation of the relation between continuous-time and discrete-time signals through sampling as the other answer outlines. Instead here I put a shorthand result.

First observe that for any integer $d$:
$$ \delta[n-d] = \frac{ \sin(\pi(n-d)) }{\pi(n-d) } = \text{sinc}(n-d)$$ is satisfied.

The righthand side is a sampled (and therefore discrete) sinc pulse whose continuous equivalent is $$\text{sinc}(x) = \frac{\sin(\pi x)}{\pi x}$$

Furthemore observe that the discrete-time sinc pulse $\sin(\pi n) /{ \pi n} $ is actually the impulse response of an ideal lowpass filter with a cutoff frequency of $\omega_c = \pi$ which is indeed an allpass filter with a frequency response of:
$$ H_{apf}(\omega) = 1 ~~~,~~~ \text{for all } ~~~|\omega| < \pi $$

As a consequence, the frequency response of a shifted sinc pulse
$ \frac{ \sin(\pi(n-d)) }{\pi(n-d) }$ will be:

$$H_s(\omega) = e^{-j \omega d } $$ for any $d$, integer or real, according to the time-shift property of the Fourier transform.

Hence it can be concluded by the inverse Fourier transform relation of the last statement that the impulse response $h_b[n]$ associated with the frequency response $$H_b(\omega) = e^{-j \omega (3/2) }$$ is
$$h_b[n] =  \frac{ \sin(\pi(n-3/2)) }{\pi(n-3/2) } $$

Note that $h_b[n]$ cannot be obtained from $ \frac{ \sin(\pi n) }{\pi n}$ by just shifting it $3/2$ samples which is not possible. Therefore you should obtain the samples of the impulse response $h_b[n]$ by evaluating its formula for each n.


Alternative to Extended Kalman Filter When Prediction Function Is Not Differentiable

I will tell you something, even if it is differntiable, use Unscented Kalman Filter for any non linear case.

This flavor of Kalman Filter, based on the Unscented Transform, is almost always superior to the Extended Kalman due to its properties.
The main reason is it is able to better predict the mean and variance (Which all Kalman Filter needs) of the distribution after the non linear function applied.

In most cases it is even simpler to calculate hence more computationally efficient.
Think of it as Poor's man Particle Filter (Though the sampling points are deterministic).


What kind of effect is this and what is it caused by?

Looks like pre-ringing. That happens whenever you apply a linear phase filter to a signal.  In your case that could be the sigma-delta converters or any type of sample rate conversion or other unknown processing that the operating system or the drivers doo.


Skeletonization of a Circle - Image Processing

The question here is: "The output of which skeletonisation"? The topological or the morphological?

It seems from what you mention on your discussion with your friend that you are referring to the morphological skeleton.

In the case of the topological skeleton, the answer is trivial. Given that the topological skeleton is the set of all foreground pixels that are at maximal distances from the boundaries of the image and the definition of the disk (all points in space that are at most at distance $R$ from the centre of the disk), the skeleton of a disk (a circle that is filled*) is the single pixel at its centre.

The morphological skeleton works slightly differently. It is based on the morphological operations of erosion followed by opening, both of which are using a structuring element.

Roughly speaking, you can think of erosion followed by opening as progressively "pealing off" boundaries from a region, until there is nothing left to peal off. (This is not entirely accurate and I will return to it further below).

Here for example is some random disc I placed on a dark background...

...and here are a few of the first iterations of morphological skeletonisation:

Where it is clear that the disc is shrinking and this looks like going towards our previous thinking regarding the topological skeleton, where eventually, all that is left is going to be a single pixel at the centre.

However, we can already see that this shape does not look so much like a circle now, it's kind of pointy and also, if we continue this process, something odd starts happening (in my case after about iteration 12 or 13):

Ultimately, what is left behind is a cross.

Where is this cross coming from?

This is a byproduct of two things:

The fact that the structuring element is discrete and its curvature at different iteration levels does not match exactly the curvature of the given disc. To understand this better, consider this extreme example that could well be showing a "disc" or a "diamond" or a square rotated 45 degrees:

$$\begin{matrix}
    0 & 0 & 1 & 0 & 0 \\
    0 & 1 & 1 & 1 & 0\\
    1 & 1 & 1 & 1 & 1\\
    0 & 1 & 1 & 1 & 0\\
    0 & 0 & 1 & 0 & 0\\
\end{matrix}$$

(In other words, the disc is not exactly a disc and it is not strictly enforcing that rule of "all pixels at distance $R$" strictly to all of its pixels. It is more like, pixels less than that distance in which case, its discrete boundary "plays" in and out of the actual continuous boundary.)

The way that the "pealing off" takes place. In fact, the morphological skeleton is the union of $n$ partial skeletons. Each partial skeleton of some image $I$ via a structuring element $B$ is obtained via:

$$ S_n(I) = (I \ominus nB) - (I \ominus nB) \circ B $$

Where $I$ is the image (here, our bigger disc), $n$ is the iteration level, $B$ is the structuring element, $\ominus$ represents erosion, $\circ$ represents dilation and $nB$ is interpreted as "$n$ applications of erosion". Also notice here that the $-$ operator has the lowest precedence. So you don't evaluate it until you have evaluated both the left-hand and right-hand expressions.

The "pealing off" is achieved by that $-$ operator in the middle but there is a little bit more into it. A partial skeleton is the difference (as in  set difference) between the "eroded" $I$mage by $B$ and the "dilation of that eroded $I$mage" by $B$. Which is basically like the difference between a shrunk image and the shrunk image enlarged. This difference is $\emptyset$ for the first few iterations but as the disc gets smaller and smaller and now the differences in the boundary start becoming more pronounced, the shrunk and shrunk-enlarged images start differing by 1 pixel (exactly at the peaks as you can see from the matrix example above). Since the skeleton is the union of the partial skeletons and each partial skeleton will have 4 pixels difference to its previous one, what "survives" throughout the iteration levels is that cross.

To understand better why I am referring to erosion and dilation as shrinking and enlarging and their difference as some "boundary", you need to understand the effect of those particular operations in detail. To do this, in addition to the excellent wikipedia articles, you can have a look at Gonzalez & Wood's "Digital Image Processing". In my second edition, chapter 9 is devoted to morphological image processing (and is totally worth a complete thorough reading) and the morphological skeletons are presented in pages 543-545.

Hope this helps.

*: The topological skeleton of a circle is a circle.


Kalman filter for tracking sinusoidal motion

Assume that the amplitude remains constant as well as the angular frequency $\omega$. The phase will be predicted using $\phi_k=\phi_{k-1} + \omega T_s$

The measurement matrix usually denoted as H and is  the Jacobian matrix or the first derivatives matrix of the model with respect to the state variables. B will be the all zeros matrix.


Can I start testing my NMF with extension on sample that have only one source?

Yes, that's a first sanity check. Basically, the idea of averaging multiple observations of a single signal buried in noise is a basic issue in signal processing. You can have several observations with multiple noise realizations and statistics, with different delays on the signal, with different data transformations (more or less linear).

If this does not work in your case, don't expect to much of NMF performance on multiple sources, and verify assumptions or test other source separation methods.


Why does the autocorrelation get its peak at zero?

Are you looking for a formal proof or the intuition behind this? In the later case: "Nothing can be more similar to a function than itself". Autocorrelation at lag $\tau$ measures the similarity between a function $f$ and the same function shifted by $\tau$. Note that if $f$ is periodic, $f$ shifted by any integer multiple of $\tau$ and $f$ coincide, so the autocorrelation has a comb shape - with peaks at the integer multiples of the period with the same height as the central peak.


Time-differentiation property of Fourier transform for $\sin(\omega_0t)$ and $\cos(\omega_0t)$

It does work. Note that the Fourier transforms you stated in your question are not the transforms of $\cos(t)$ and $\sin(t)$, but of $\cos(\omega_0t)$ and $\sin(\omega_0t)$.

Let's Assume we know the transform of $\sin(\omega_0t)$:

$$\mathcal{F}\{\sin(\omega_0t)\}=\frac{\pi}{j}[\delta(\omega-\omega_0)-\delta(\omega+\omega_0)]\tag{1}$$

Now we take the derivative of $\sin(\omega_0t)$:

$$[\sin(\omega_0t)]'=\omega_0\cos(\omega_0t)\tag{2}$$

Using $(1)$ and the differentiation property of the Fourier transform, the Fourier transform of $(2)$ must be

$$\mathcal{F}\{\omega_0\cos(\omega_0t)\}=j\omega\cdot \frac{\pi}{j}[\delta(\omega-\omega_0)-\delta(\omega+\omega_0)]\tag{3}$$

Noting that for any function $f(\omega)$ that is continuous at $\omega=\omega_0$ we have

$$f(\omega)\delta(\omega-\omega_0)=f(\omega_0)\delta(\omega-\omega_0)\tag{4}$$

and, consequently, we can rewrite $(3)$ as

$$\begin{align}\mathcal{F}\{\omega_0\cos(\omega_0t)\}&= \pi[\omega_0\delta(\omega-\omega_0)-(-\omega_0)\delta(\omega+\omega_0)]\\&=\omega_0\pi[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)]\end{align}\tag{5}$$

as expected.


CDMA scheme that preserves constant envelope (NRZ/ASK shape)

Yes you are correct that the reception of multiple CDMA signals from different transmitters tends toward a complex Gaussian distribution given the central limit theorem at play - which is not very different from the thermal noise that a spread spectrum signal may be buried in at reception. Yet this does not preclude us from making a hard decision on each sample at the receiver! So there is no need for a constant envelope CDMA in order to use a simple slicer.

Since CDMA signals are buried in the "noise" of other signals and natural noise, if you take multiple samples and average you can get higher effective precision depending on the number of times you over sample within any symbol interval. Oversampling and averaging in the presence of noise is called "dithering"; as long as there is zero-mean noise present the average will converge to higher precision. (Simple example: consider a 1 bit converter that puts out 0 or 1 based on the input being higher or lower than 0.5. If truth was .25, and if we had a uniform noise distribution on each sample that exceeded +/-0.5, then over many samples we will decide 3 zeros for every 1 zero and the average with be 0.25.)

If you do this simplistically by just increasing your sampling rate, you end up getting an extra $0.5$ bits of precision for every doubling of the sampling rate. See this post for further details: What are advantages of having higher sampling rate of a signal?

Sigma Delta converters do this with noise shaping for even higher gains in precision vs the sampling rate.

Given the reception of CDMA is to correlate by accumulating the results of each sample over long codes before deciding what symbol was actually sent (see How CDMA receiver extract it's corresponding data from the receiving modulated & superposition-ed signal?) , it is not even necessary to oversample to achieve the same effect of averaging by dithering as described above- the noise gets averaged out in the correlation.

This is demonstrated below as applicable to both a standard bi-phase (BPSK) communication and BPSK with CDMA with the descriptions above. In the plot we see a barker code sequence [1 0 1 1 0 1 1 1 0 0 0] which has the desirable property as a PN sequence to not be correlated to shifted versions of itself, which was detailed at the link given above. As transmitted it is either sent as a $+1$ or a $-1$. As received it is corrupted by noise and the histogram of all received signals given a $+1$ was transmitted (blue) or a $-1$ was transmitted (red) is shown on the left. For the samples in the histogram that go beyond the zero-crossing, an error would result when we make a decision on that sample. In standard BPSK communication we are indeed concerned about the errors in every bit, to the amount of that histogram crossing the zero-crossing is a big deal and results in the waterfall curves we get relating BER to SNR (as noise in the SNR is the standard deviation of that histogram related to S which is the average).

HOWEVER notably in CDMA any individual error is of less concern since we are no longer concerned with the individual decisions (hence we will call the individual values in the sequence chips and not bits) but the overall average of our decisions over the complete symbol interval.

Consider the two plots below, which is the same me as our plot above showing a BPSK sequence in high and low SNR conditions. The upper plot with what would result at a higher SNR and a relatively lower error rate on each chip. The bottom plot shows the typical CDMA application where we we receive the signal in the presence of very high noise, including the signals from other CDMA transmitters using other codes (so to us they appear very much as AWGN noise). In this case, the tail showing the distribution of samples that would result in an error on each chip is quite large. However to be clear in CDMA we are capturing lots of samples before making a final decision as to which overall code was transmitted--an individual chip error is not a significant consequence. So in the lower plot as long as we are averaging a very large number of samples before we make our final decision (In GPS there are 1023 chips in each sequence, and we also go through 20 sequences before making a data bit decision, so 20460 samples if we only used one sample per chip) we will still demodulate with a very low error rate.

If that is still confusing consider that we undo the code sequence by multiplying it with the code in the receiver, so if a data symbol of one is desired, the transmitter will send the sequence as [1 -1 1 1 -1 1 1 1 -1 -1 -1], and if a data symbol of of zero is desired the transmitter with send the inverse of the code as [-1 1 -1 -1 1 -1 -1 -1 1 1 1]. Thus we convert the code sequence to all one's for data symbol "one" and all minus one's for data symbol "zero". In the presence of AWGN noise only a simple repetition code would have the same effect.

However to note: A typical bi-phase CDMA signal (such as GPS) given the typical phase, time and frequency offset of the receiver will be equally distributed in AM and PM components. Slicing is equivalent to hard-limiting which removes all AM. You may conclude from this that there would be a 3 dB loss of signal; actual results range from 1.6 to 2.2 dB depending on the excess bandwidth of the receiver prior to slicing (more bandwidth less noise penalty). There is a nice chart showing this in James Spilker's GPS books "Global Position Systems Theory and Practice Volume I and II".

This is illustrated further with the graphic below depicting operation of CDMA in the presence of a non-coherent single-tone jammer where in fact there would be a 3 dB loss in SNR under this condition due to hard-limiting the signal. This helps illustrate the important distinction on the effects of hard limiting a CDMA signal that we would expect to receive with a negative SNR, compared to a traditional BPSK signal that we would expect to receive with a positive SNR.

For high SNR conditions with the signal much larger than the jammer, the diagram on the left shows the signal as received with some arbitrary phase offset (and would also have a carrier offset, so this graphic would be rotating at that offset rate). The fixed phasor represents our signal relative to the jammer that is shown in the smaller phasor rotating with its relative offset in frequency (non-coherent). Here the noise due to the jammer is composed of equal AM and PM components, and with hard limiting we remove all the AM without changing our signal level resulting in a +3 dB improvement in SNR.

For low SNR conditions, the situation is reversed in that our signal is much smaller than our noise which is the typcial case for CDMA operation. So here our signal prior to hard-limiting is the AM and PM components on the single tone jammer. With hard limiting we reduce our signal level 3 dB further relative to the jammer, thus is a loss of 3 dB in SNR.

The other issue with hard-limiting and slicing is vulnerability to coherent jammers. If your receiver is slicing and an intentional jammer is able to send a single tone coherently in phase then then received signal will be an AM only component on top of the jamming signal and will be completely removed.  Jamming and interference rejection in general is a high motivation to use higher levels of quantization for CDMA signals; there is little benefit of going beyond 3-4 bits otherwise. So for AJ considerations it is a really bad idea to hard-limit or equivalently slice.


What is the $\mathcal Z$-transform of Bessel function $J_0(\alpha n)$ sequence

The Taylor expansion for the Bessel function of the first kind and 0th order is
$$
J_{0}\left ( x \right )= \sum_{m=0}^{\infty }\frac{(-1)^{m}}{(m!)^{^{2}}}\left ( \frac{1}{2}x \right )^{2m}
$$

(see http://en.wikipedia.org/wiki/Bessel_function)

So you can basically approximate this as the Z-transform of a polynomial.


Compressive Sensing - Sparse in frequency example

I think I understand your confusion. In your first example, it would have been more realistic to take

y = x(q)

for example.

A possible explanation for the approach you see here is that authors of such examples would like to ensure they obtain the exact sparse solution they generate the example from to compare against. Therefore they generate the measurements directly from it. It is not "cheating" per se in a mathematical sense, but you can say that it does not reflect the practical reality of taking measurements in the time domain.

What could happen when taking the measurements in the time domain is that if you generate sinusoids with arbitrary frequencies, the signal will appear less sparse in the frequency domain when the frequencies do not match the frequencies of the dictionary. If for example you generate a sinusoid at 1.5 Hz and your DFT has atoms corresponding to the frequencies 1 Hz, 2 Hz, 3 Hz, etc., your signal will not be 1-sparse in this dictionary but will have to use multiple atoms at adjacent frequencies to form the sinusoid at 1.5 Hz. A simple way to avoid this is to generate a signal of the prescribed sparsity directly in the frequency domain and then convert it to the time domain as part of taking the measurements.

Technically, it does not involve already knowing the frequency spectrum of the signal, because X is not known in the reconstruction and the sufficient frequency information in the classic sampling-theoretical sense is lost in the multiplication

y = A*X

because A is a wide matrix.


My first question is, why should we allow the signals have different phases?

Counterquestion: Why should they not be allowed? Enforcing continuous phase can make systems more complex, without any benefit.

I cant see a reason for this and as far as I know, different phases result in phase discontinuity of the transmit signal.

Depending on how you do that transition, that's not necessarily the case: you can pulse-shape $\theta_i$, or you can filter $x$.
It's less often a problem than old textbook make it out to be, to be honest.

(2) And the second question is, how this difference in phase is related to coherent and non-coherent receiver?

Not at all. What the receiver does says nothing about what the receiver does!


BlackMath hit the nail on its proverbial head:

  Using control channels from the base station.

Typically, the whole circus works like this:

Before a user equipment can register at a base station, it has to look for that cell's "beacon". Whenever it sees that, it waits a time $T_\text{wait}$ after the beacon, at which new devices are allowed to say "hey, I'm here". Because these devices can't know of each other yet, that phase isn't really TDD per se (it's more like ALOHA, with the random element being delay do to distance); especially, because these devices are at different distances from the base station (which they can't know already), it's not time-slotted. For that reason, the "hello" phase is not as tightly timed as regular TDD slots. That means it "wastes" a lot of uplink time, and for that reason, the network operator might define these phases to happen relatively rarely.

Assuming the new user equipment successfully registered within that "hello" phase (i.e. didn't collide with another new device), the base station now knows that there's a new device, and this devices cumulative delay $T_\text{beacon} + T_\text{distance} + T_\text{wait} + T_\text{distance}$. Since to the base station, $T_\text{beacon}$ and $T_\text{wait}$ are known, it can easily deduct $T_\text{distance}$. It then uses the downlink to tell the new device that, "hey, you need to be $T_\text{timing advance}$ earlier than you are, to fit into our TDD scheme. Your timeslots to send data are 1, 3, 22 and 128. Welcome to my cell!".

The user equipment adjusts its internal TDD timekeeping using that timing advance. During operation, the cell tower continuously informs the user equipment whether it needs to do adjustments of that timing advance, which would, for example, be the case if the user moved away or towards the cell tower.


okay, using the snipping tool, i tried to put in just the 2nd-order IIR lattice in the following drawing.  $x[n]$ is being input in the upper left.  the numbering for the $v_*$ coefficients is horrible, but i just cut this drawing out from somewhere else.

now the transfer function i get for this is

$$ H(z) = \frac{Y(z)}{X(z)} = \frac{v_4 + v_3(-k_1 + z^{-1}) + v_2\big(-k_2+k_1(k_2-1)z^{-1} + z^{-2}\big)}{1 + k_1(k_2-1) z^{-1} - k_2 z^{-2}}$$

and if you compare that to the standard direct form

$$ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}} $$

you should have equivalence when

$$\begin{align}
k_2 &= -a_2 \\
k_1 &= -\frac{a_1}{a_2+1} \\
\\
v_2 &= b_2 \\
v_3 &= b_1 - a_1 b_2 \\
v_4 &= b_0 - \frac{a_1}{a_2+1} b_1 + \left( \frac{a_1^2}{a_2+1} - a_2 \right)b_2 \\
\end{align}$$

i would appreciate it if someone else would check my work.  i could well have fucked up a sign or something.

so, wherever you get your five Direct Form $a_k$ and $b_k$ coefficients (e.g., the Cookbook, but remember to divide everything by $a_0$), you can convert them to the five Lattice Form coefficients with the above equations (assuming i didn't mess it up).


Image Contrast Enhancement based Sub-histogram Equalization Technique without Over-equalization Noise

The paper is not really well written, so it's not surprising that you're having a hard time. These are a few hints to guide you:

the authors try to do some mode detection on the image brightness histogram in order to segment it between common values in the image (= high density in their terminology) and rare values in the image (= low density);
in order to segment the histogram, they fix a window width (i.e., a number of bin of the histogram), apply this window over all possible positions in the histogram and sum the populations of the bins inside the window ($Eq. (7)$). This gives the density of a histogram bin (or equivalently of a brightness level);
their segmentation process consists then in going from the peaks of the histogram and merging neighbouring bins into it as long as the difference between the border bin and its neighbour (= the momentary change) is less than the average of the momentary changes in the region.


How to evaluate a background subtraction algorithm given ground truths?

One way would be to create a simple signal-to-noise estimate, using the pixel intensity and the ground truth image as a mask. Here's an example of what I mean in Python:

import numpy as np

# Sample data
# Indices [0, 1, 2, 6] are background, [3, 4, 5] are object
original  = np.array([1, 2, 2, 5, 6, 8, 2])
truth     = np.array([0, 0, 0, 1, 1, 1, 0])
back_mask = 1 - truth
result    = np.array([0, 1, 0, 4, 5, 7, 1]) # Result of background subtraction

def non_masked_mean(input_array, mask):
    """mean of non-masked elements in the array"""
    return np.ma.masked_array(input_array, mask).mean()

def snr(input_array, back_mask):
    """ratio of non-masked to masked mean intensities"""
    with np.errstate(all='ignore'):
        snr = non_masked_mean(input_array, back_mask) / non_masked_mean(input_array, (1-back_mask))
    return 0 if np.isnan(snr) else snr

snr_before = snr(original, back_mask)  # (mean of objects / mean of background), before
snr_after  = snr(result, back_mask)    # (mean of objects / mean of background), after
snr_ratio  = snr_after / snr_before

if snr_after == np.inf:
    print "Perfect SNR in result! Either great or suspect..."
elif snr_after == 0:
    print "Image was flattened :("
else:
    print "SNR changed by a multiple of %.2f" % (snr_ratio)

The problem with this is that you get a perfect score by favouring a heavy-handed background subtraction. Even if only 1 pixel of the object survives this SNR will be infinite. Another approach would be to consider the background and the signal separately, and weight the results of the two. Then you can decide if it's more important to flatten the background, or preserve the original.

def non_masked_mean_ratio(input_array, result_array, mask):
    """ratio of two array means, with the same mask"""
    return non_masked_mean(result_array, mask) / non_masked_mean(input_array, mask)

signal_ratio     = non_masked_mean_ratio(original, result, back_mask)     # (mean of objects after) / (mean of objects before)
background_ratio = non_masked_mean_ratio(original, result, truth)         # (mean of background after) / (mean of background before)
weight           = 0.5   # Higher weight favours preservation of objects

score = (weight*signal_ratio) + ((1-weight) * (1-background_ratio))

This gives the score a nice expected range of about 0-1, where 1 is the ideal result. It's not perfect though, as you can boost it by increasing the intensity of the objects.

If you're expecting a binary result, you can use a similar strategy except using the number of non-zero pixels instead of the mean:

def nonzero_count_ratio(result_array, ideal_array):
    """ratio of actual number of non-zero pixels vs ideal non-zero pixels"""
    return float(np.count_nonzero(result_array*ideal_array)) / np.count_nonzero(ideal_array)

true_obj_ratio   = nonzero_count_ratio(result, truth)     # Number of correct object pixels vs number of potential correct
false_back_ratio = nonzero_count_ratio(result, back_mask) # Number of false background pixels vs number of potential false
weight           = 0.5

score = (weight*true_obj_ratio) + ((1-weight) * (1-false_back_ratio))


Note that the second half of the period of the sinusoid is simply equal to the first half times $-1$. Since the same is true for the square wave, the multiplication of the two is the same in both intervals $[0,T/2]$ and $[T/2,T]$. That's why the integral over the whole period equals twice the integral over the first half of the period.


Get the transfer function when using the "designfilt" command in MATLAB

The filter is a discrete FIR filter, so the transfer function would be given as:
$$H(z) = .0287 + 0.1430z^{-1} + 0.3283z^{-2} + 0.3283z^{-3}+0.1430z^{-4}+0.0287z^{-5}$$
With the frequency response given by using the unit circle for the complex variable $z$, as in $z = e^{j\omega}$ for $\omega \in [0, 2\pi)$, with the sampling rate normalized to $\omega = 2\pi$.
You can use the freqz command in Matlab to plot or return the frequency response:
freqz([0.0287 0.1430 0.3283 0.3283 0.1430 0.0287])

Which creates a plot similar to what I have done using Octave below for these coefficients:

The frequency axis will be normalized to the radian angular sampling frequency, but you can refer to Matlab's documentation for other parameters that can be passed in to freqz to scale the frequency otherwise.


DSD bits are 0 and 1 indeed, but they actually represent normalized outputs of -1 and 1, respectively, since one-bit samples can only be mid-raiser-quantized. This means DSD samples cannot represent zero. One can generate DSD silence by alternating 0-s and 1-s (your second option), or by a 1-bit delta-sigma modulation of a zero PCM input, with or without dithering (your third option).
Multi-bit modulators can (and normally use) mid-tread quantization, so they can represent zero. Thus, multibit silence is either a sequence of zeroes or the output of a multi-bit delta-sigma modulator presented with a zero input.


Seamless looping of a signal without pops

But is this good practice?

Yes, that's a perfectly reasonable approach.

Will this will work for all types of signals I might have ?

There are a LOT of signals out there and there is always an outlier. In your case, I think the most "vulnerable" would be low frequency sine waves. A 40 Hz sine wave has a period of 25 ms and looping one at an inconvenient point without audible artifact is tricky. So this all boils down to: how likely are these signal for your application, what level audibility is acceptable and what are the constraints around resources,  looping points & length, etc.

Is there any other way this could be done in a more "proper" manner?

If you have some leeway with the exact loop length and loop points, try to put the loop points at the nearest zero crossing. Ideally both zero crossing go in the same direction.
You can do an energy cross fade instead of a linear cross fade. This will minimize the volume dip in the middle of the cross fade, but at 15ms it's not going to make much of an audible difference.


get convex hull of isolated shapes with openCV/C++

Try dilation, a technique used to enlargen light spots over darker spots in an image.

Start by dilating the image such that many of the enlargened white blobs will overlap and form a complete connected contour, as such:

Now, when OpenCV finds contours from this large figure to make convex hulls, it will find a complete circle rather than loose seperate blobs.

From here, you can find the largest contour, which is the (outer) largest light-purple contour.

Here is the code for this (which currently computes the contours, but not the largest area - this is just to demonstrate the dilation):

#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
#include <stdio.h>
#include <stdlib.h>

using namespace cv;
using namespace std;

int main( int argc, char** argv )
{
 Mat src; Mat src_gray;
 src = imread( "convex.jpg", 1 );
 //resize(src, src, Size(640,480), 0, 0, INTER_CUBIC);
 cvtColor( src, src_gray, CV_BGR2GRAY );

 namedWindow( "Source", CV_WINDOW_AUTOSIZE );

 //Dilate
 int erosion_size =10;
  Mat element = getStructuringElement(cv::MORPH_ELLIPSE,
        cv::Size(2 * erosion_size + 1, 2 * erosion_size + 1),
        cv::Point(erosion_size, erosion_size) );

 dilate(src_gray,src_gray,element);

 imshow( "Source", src_gray );

 // Convex Hull implementation
 Mat src_copy = src_gray.clone();
 Mat threshold_output;
 vector<vector<Point> > contours;
 vector<Vec4i> hierarchy;

 // Find contours
 threshold( src_gray, threshold_output, 200, 255, THRESH_BINARY );
 findContours( threshold_output, contours, hierarchy, CV_RETR_TREE, CV_CHAIN_APPROX_SIMPLE, Point(0, 0) );

 // Find the convex hull object for each contour
 vector<vector<Point> >hull( contours.size() );
 for( int i = 0; i < contours.size(); i++ )
 {  convexHull( Mat(contours[i]), hull[i], false ); }

 // Draw contours + hull results
 RNG rng;
 Mat drawing = Mat::zeros( threshold_output.size(), CV_8UC3 );
 for( int i = 0; i< contours.size(); i++ )
 {
  Scalar color = Scalar( rng.uniform(0, 255), rng.uniform(0,255), rng.uniform(0,255) );
  //drawContours( drawing, contours, i, color, 1, 8, vector<Vec4i>(), 0, Point() );
  drawContours( drawing, hull, i, color, 1, 8, vector<Vec4i>(), 0, Point() );
 }

 // Show in a window
 namedWindow( "Hull demo", CV_WINDOW_AUTOSIZE );
 imshow( "Hull demo", drawing );

 waitKey(0);
 return(0);
}


Understanding 2D FFT result of an image having a pattern

I found that FFT output depends on aspect ratio of the input image. Fourier pattern is perpendicular to the input image pattern only if the input image has same width and height.


How to show that y[n] = x[n] * h[n] turns into the Y(z) = X(z).H(z)?

The correct sequence $y[n]$ is the one you obtained via the $\mathcal{Z}$-transform:

$$y[n]=\frac{1}{a-b}\left[a^{n+1}-b^{n+1}\right]u[n]\tag{1}$$

I'm not sure where that $n+1$ term is coming from in your time-domain derivation, so I'll write out the convolution sum:

$$y[n]=\sum_{k=-\infty}^{\infty}a^ku[k]b^{n-k}u[n-k]=u[n]\sum_{k=0}^{n}a^kb^{n-k}=u[n]b^n\sum_{k=0}^{n}\left(\frac{a}{b}\right)^k=\\=u[n]b^n\frac{1-\left(\frac{a}{b}\right)^{n+1}}{1-\frac{a}{b}}=u[n]b^{n+1}\frac{1-\left(\frac{a}{b}\right)^{n+1}}{b-a}=u[n]\frac{b^{n+1}-a^{n+1}}{b-a}$$

which equals (1), as it should.


Measuring frequency response range of piezoelectric disc

I am afraid that it is rather impossible without a proper hardware. Sweep sine is ok as a general method, but you would either need:

Reference transducer with known (preferably) linear frequency response. Then you can find the difference between those two.
Signal actuator with known frequency response, and then you can also find the difference between recorded signal and what it should be.

In any other case you will only measure the convolution of actuator response and response of piezoelectric element.

You should also understand that this disc is meant to be used as a contact device, not for recording of sound waves. There is a huge difference between those two signals recorded.

Lastly you can try to search for frequency characteristic of known model (presumably measured by someone). I did similar work back in time, but we used different type.


What's labeled $(1)$ in your question is a special case, the flat channel. It can be represented as a single coefficient.

In general, channels aren't flat, and we then need to apply $(2)$ instead. That's no different from acoustics.

So, your claim that $(1)$ generally applies to SISO channels is plain wrong.

However, when naming something "SISO", one usually means "in contrast to MIMO"; and in the context of MIMO, we basically need our channels to be flat, so that the MIMO channels can be represented as a $m\times n$ matrix of single coefficients between $m$ and $n$ transmitting/receiving antennas, respectively.¹

But, you say,

  if the individual channel is like $(2)$, how can we model it like $(1)$?

We can't. But we can apply multi-carrier techniques (mostly, OFDM), to first divide a wideband channel with a convolutive channel impulse response into many narrow channels, which individually are flat.

¹ In theory, one can do MIMO on multipath / frequency-selective channels instead of first dividing them in many flat subchannels, but as far as I know that's not done commercially. OFDM implementations are mature and cheap, and it's much easier to come up with working math that allows us to decompose the superpositioning MIMO channel when it's flat. In fact, there's no standard textbook I'm aware of that even does MIMO on frequency-selective channels.


What is an intuitive understanding of the difference between the minimum required bandwidth and actual bandwidth?

We must be careful in how we define bandwidth. To maintain truly rectangular pulses in the time domain, infinite bandwidth is required, as in the frequency domain the spectrum would be a Sinc function that extends to infinity (as the Fourier Transform of the rectangular pulse). The requirement for infinite bandwidth is intuitively consistent with a rectangular pulse that can change instantaneously from a low level to a high level (infinite bandwidth is required to transition that fast!). Note for now that the Sinc function has nulls at 1/T and integer multiples of that, where T is the width of the rectangular pulse.

As we constrain the bandwidth in frequency, the edges of the signal will likewise slow in transition from low to high. A good rule of thumb related to this (specific to first order systems but generally a good estimate) is:

$$t_r = 0.35/BW$$

Where $t_r$ is the 10% to 90% rise time in seconds, and BW is the bandwidth in Hz.

Ultimately for a digital signal the fastest bandwidth is required when we send the sequence 1 0 1 0 1 0 1 0 ....   We note that for rectangular pulses this would be a 50% duty cycle square wave which according to the Fourier Series Expansion is decomposed in frequency to the odd harmonics of the fundamental sinusoidal frequency that fits in the 1 0 1 0 pattern.  This is still the Sinc function we described above, just when the time domain signal repeats the only frequency components that can exist are at integer multiples of the repetition rate, and therefore those frequencies become samples of the Sinc function.  As described above, the nulls would exist at 1/T, where T in this case is $1/R_b$ which is also the location of the 2nd harmonic that we already reasoned is nulled (only odd harmonics exist for the 50% square wave). The first harmonic is at $R_b/(2)$, and if we were to represent it as a double sideband signal (positive and negative frequencies) then the minimum bandwidth required to contain both would be $R_b$. This would result in a sinusoidal waveform (sine or cosine) but would be sufficient to represent our 1 0 1 0 pattern without data loss (we recover with a simple threshold at y=0).

However, we can also represent the same sequence as a single sided signal (positive or negative frequencies only) which would then only require half that bandwidth, which is how you came to $R_b/2$. This would be consistent with a use of Single Sided Modulation.  (Note we don't require any bandwidth if the pattern was always 1 0 1 0 as it would be represented by a single frequency, but this is the case of the highest frequency needed, and we would also need all the lower frequencies to represent any other patterns).

Depending on your background this may be a lot to swallow, and if so I also point you to the links below that also help explain these fundamental concepts:

How do I modulate complex time domain signal to a carrier frequency?

Frequency shifting of a quadrature mixed signal


Note that the inverse of an FIR system is IIR, and the same is true for the inverse of an IIR system, unless it is an all-pole system, the inverse of which would be FIR. So in most cases the ideal equalizer should have an infinitely long impulse response in order to perfectly invert the channel. In practice almost all adaptive equalizers are FIR filters which can only partly compensate for the channel response. This is why I think that there should be no reason for your equalizer to perform worse for IIR channels than for FIR channels. You could post some more information about the channel responses and the convergence of your equalizer to make it easier to see what's going on.

As for the zero forcing criterion, such criteria for designing an equalizer only make sense for noisy channels, because otherwise there is no trade-off between equalizing the channel response and noise enhancement. Zero-forcing and MSE criteria result in the same equalizer in the absence of noise. If you have no noise, the only problem is the (hopefully linear) distortion of the channel and naturally this is the only thing the equalizer would compensate for. However, with an FIR equalizer exact zero-forcing is usually impossible as discussed above.


for overlap-save, (sometimes called "overlap-scrap") when you (circularly) convolve h with x1. you scrap the first 3 output values of y1 and keep the last 5 values.  when you (circularly) convolve h with x2, again you scrap the first 3 values of y2 and keep the last 5 values.

with

x[n]={7,6,4,5,2,4,5,2,3}
h[n]={1,2,3,1}

it should be

x1[n]={0,0,0,7,6,4,5,2}
x2[n]={4,5,2,4,5,2,3,0}

remember that in convolution, h[n] gets turned around to be h[-n] (h[0] remains in the same position but the other values are extended to the left).  only when h[0] is aligned with x1[3] will you get a legitimate value for y1[3].  and notice, because you (necessarily) delayed x[n] by 3 samples with that initial zero-padding, the output will also be delayed by the same amount.


You are absolutely right. With $f=F/F_s$ you obtain the equality

$$\int_{\frac12}^{\frac12}X(f)e^{j2\pi nf}df=\frac{1}{F_s}\int_{-F_s/2}^{F_s/2}X\left(\frac{F}{F_s}\right)e^{j2\pi nF/F_s}dF\tag{1}$$

which is also what I find in my third edition of the book (the right-hand side of $(1)$ equals the left-hand side of Eq. $(4.2.81)$ in my edition).


The solution turned to be the right number of overlapping input buffers.
To get the source into fft and then back we need to apply windowing on the sources to reduce frequency leaks. Reconstructing such signal without any modification works just fine. However as @applesoup poined out certain operations in the frequency domain introduce their own discontinuities at the edge of output buffers so the output buffers must be windowed as well before combining into the output signal which in turn introduces inconsistency to the amplitude of the signal as in effect we need to apply same window twice. The solution to this problem is to double the number of buffers and change the overlap to 25% (of looking from the other side, to 75%) so each sample is reconstructed from 4 windowed output buffers which compensates the double windowing of each individual buffer.


best technique to remove noise from bio signal data

As a first appraoch I would suggest to use a mean filter, which will calculate the average of the last n values.

$$y[n] = 1/N \sum_{k=0}^{N-1} x[n-k]$$

If you have access to the full signal (offline processing) you can perform forward and backward filtering which will avoid the delay introduced by a single filter. (Check matlab's filtfilt function for example)


Non-symmetrical DFT of real valued signal

Your choosen example frequency $\omega_0 = \frac{5 \pi} {32}$ touched a numerically sensitive case, in which the imaginary part of the theoretical DFT which should be zero is not observed so, on the numerical computation. Since the imaginary part is computed as non-zero, the phase which should be zero according to $\phi[k] = \tan^{-1}( \frac{\Im\{X_1[k]\}}{\Re\{X_1[k]\}} )$ whenever $\Im\{X_1[k]\}$ is zero, turns out to be some erratic function of $k$ (Note also that real part of $X_1[k]$ is also not computed to be exactly zero for those $k$ for which it should be, which also contributes to the erratic phase computation behaviour)

Consider the discrete time signal $x_1[n] = \cos(\frac{5\pi n}{32})$ and its 64-point DFT $X_1[k]$ (for ineteger $k \in [0:63]$) which is computed as:

$$ X_1[k] = \sum_{n=0}^{63} \cos(\frac{5\pi n}{32}) e^{-j \frac{2\pi}{64} nk} $$

expand the cosine into complex exponentials:
$$ X_1[k] = 0.5 \sum_{n=0}^{63} ({e^{j \frac{5\pi}{32} n} + e^{-j \frac{5\pi}{32} n}}) e^{-j \frac{2\pi}{64} nk} $$

adjust the cosine based complex exponential arguments:
$$ X_1[k] = 0.5 \sum_{n=0}^{63} ({e^{j \frac{10\pi}{64} n} + e^{-j \frac{10\pi}{64} n}}) e^{-j \frac{2\pi}{64} nk} $$

separate the sums:
$$ X_1[k] = 0.5 \sum_{n=0}^{63} {e^{j \frac{10\pi}{64} n} e^{-j \frac{2\pi}{64} nk} + 0.5 \sum_{n=0}^{63}e^{-j \frac{10\pi}{64} n}} e^{-j \frac{2\pi}{64} nk} $$

Now replacing numerical value of $N = 64$ with $N$ we shall see the summations:
$$ X_1[k] = 0.5 \sum_{n=0}^{N-1} {e^{j \frac{10\pi}{N} n} e^{-j \frac{2\pi}{N} nk} + 0.5 \sum_{n=0}^{N-1} e^{-j \frac{10\pi}{N} n}} e^{-j \frac{2\pi}{N} nk} $$

Merge the complex exponentials:
$$ X_1[k] = 0.5 \sum_{n=0}^{N-1} {e^{j \frac{2\pi n (5 -k)}{N}}+ 0.5 \sum_{n=0}^{N-1}e^{-j \frac{2\pi n (5+k)}{N}}} $$

Now the above sums add up to $N/2$ whenever $k=5$ or $k=-5$ ( this is also for $k = -5+ 64 = 59$) and they become zero for any other integer $k$ in the range $k \in [0:63]$

As you can now see, $X_1[k]$ is a real and even symmetric function of k. And its imaginary part is zero, however in your loop computation this
last condition is not realized due to numerical issues and therefore you observe erratic phase computation behaviour.


Applying a 2D band pass filter (Butterworth ?) to an array with NaNs

filtfilt() is a technique to achieve zero-phase filtering by applying the same filter twice to the data; with the output of the first stage reversed and filtered again in the second stage. Zero phase filtering is a desired property in image processing.
NaN means "not a number" and indicates those indeterminate conditions like $0/0$, $\infty/\infty$, $\infty - \infty$. NaNs cannot be processed arithmetically, hence you should remove all the NaN samples with some suitable values before further processing.


Frequency resolution and timestep in DFT

The "timestep" $dt$ (using your nomenclature) between DFT windows can really be anything you choose. There are a few different cases to consider:

$dt < N$: In this case, there is overlap between successive windows of samples that you input to the DFT. This can be used if you would like to get answers more frequently than every $N$ samples. It's quite common to see something like 50% overlap between windows being used.
$dt = N$: In this case, there is no overlap; you just take successive blocks of $N$ samples and use those as inputs to the DFT. This would give you the minimum processing load while not skipping any input data.
$dt > N$: As you noted, in this case, you're skipping data between windows (sometimes called "underlap"). This can reduce your processing load further, but since you're skipping data, you run the risk of missing the analysis of some feature in that portion of data that was discarded.

The frequency spacing of the DFT is related to the sample rate and the DFT length. Specifically:

$$
df = \frac{f_s}{N}
$$

where $f_s$ is the sample rate of the input data. $df$ in the equation above corresponds to the frequency spacing between each output (i.e. with $f_s$ measured in Hz, the outputs will have a spacing of $df$ Hz).


A tapped-delay-line is really a methodology or implementation. An FIR filter is a more general description of a purely feedforward digital network.

While your impulse response is finite (no feedback components) you can consider the tapped-delay-line filter to be a particular FIR filter whose coefficients are non-zero at tap positions and zero otherwise.

Should you introduce feedback components, I think you could still describe the system as a tapped-delay line filter, but it would cease to be an FIR filter at this point and become an IIR filter.


Transfer function and difference equations: why does $H(z)$ numerator polynomial not correspond to $Y(z)$?

It is right that $H(z)=\frac{Y(z)}{X(z)}$ and yes, the output is in the numerator. But the numerator of one side is multiplied by the denomarator of the other side to get an equality:

$$H(z)=\frac{b_Mz^M+\cdots+b_0}{a_Nz^N+\cdots+a_0}=\frac{Y(z)}{X(z)}\tag{0}$$
$$\hspace{-1.5cm}\Rightarrow\left(b_Mz^M+b_{M-1}z^{M-1}+\cdots+b_0 \right)X(z)=\left(a_Nz^N+a_{N-1}z^{N-1}+\cdots+a_0\right)Y(z)\tag{1}$$
Taking the inverse $z$ transform from $(1)$ yields the differnce equation of the system:
$$b_Mx[n+M]+b_{M-1}x[n+M-1]+\cdots+b_0=a_Ny[n+N]+\cdots+a_0\tag{2}$$
So indeed, I explained it in the reverse order. Everything begins from the difference equation of a system in $(2)$, then the $z$-transom is applied on it:
$$\hspace{-1cm}b_Mz^MX(z)+b_{M-1}z^{M-1}X(z)+\cdots+b_0 X(z)=a_Nz^NY(z)+a_{N-1}z^{N-1}Y(z)+\cdots+a_0Y(z)\tag{3}$$ and then factoring $X(z)$ on the left and $Y(z)$ on the RHS, we will get to $(1)$. From that, $(0)$ is resulted.

If we [wrongly] had
$$\frac{b_Mz^M+\cdots+b_0}{a_Nz^N+\cdots+a_0}=\frac{X(z)}{Y(z)}\tag{4}$$ then after multiplication we would  get
$$\left(b_Mz^M+b_{M-1}z^{M-1}+\cdots+b_0 \right)Y(z)=\left(a_Nz^N+a_{N-1}z^{N-1}+\cdots+a_0\right)X(z)\tag{5}$$
$$\hspace{-1.5cm}b_Mz^MY(z)+b_{M-1}z^{M-1}Y(z)+\cdots+b_0 Y(z)=a_Nz^NX(z)+a_{N-1}z^{N-1}X(z)+\cdots+a_0X(z)\tag{6}$$
and taking the inverse $z$-transform we would have
 $$\hspace{-2cm}b_My[n+M]+b_{M-1}y[n+M-1]+\cdots+b_0=a_Nx[n+N]+\cdots+a_0\tag{7}$$
which is obviously not the correct difference equation (i.e. $(2)$ that we started from).


Zero padding and 2D Fourier transforms: how does zero-padding affect phase?

Symmetric zero-padding (in the center of an image around the N/2,N/2 sample) does not affect the FFT phase result.  Or after an 2d fftshift before the 2DFFT, symmetric zero-padding around the edges (circularly around the 0,0 sample) does not add phase shift.

An FFT phase measure the even to odd ratio around 0,0.  This won't change with any padding that maintains the identical even and odd decomposition of the input.  An fftshift moves the phase reference to the middle of a 2D image, which may make more sense with reference to an odd and even function decomposition of the original input.


Get position of different frequencies in sound signal

I don't understand the interferometer setup but have a look at:

Eric Jacobsen & Richard Lyons. DSP Tips & Tricks: The Sliding DFT. IEEE Signal Processing Magazine, March 2003, pp 74–80.

Eric Jacobsen & Richard Lyons. DSP Tips & Tricks: An Update to the Sliding DFT. IEEE Signal Processing Magazine, January 2004, pp 110–111.

You can use the sliding DFT (discrete Fourier transform) to calculate a windowed DFT starting at the time of each sample, not for all frequencies but for just the frequency you are interested about. See the subsections about generalized cosine windows in Wikipedia for the kind of windows that work with sliding DFT. Perhaps you can then use the magnitude and phase information from sliding DFT for your task.


Linearity and shift-invariance of 2-D system on lattice

Well, actually you only need to insert the system $H$ into the definition of linearity and check, if it holds:

$$\begin{align}H\{af_1[x]+bf_2[x]\}&=\sum_y\phi[y](af_1[x-y]+bf_2[x-y])\\
&= a\sum_y\phi[y]f_1[x-y]+b\sum_y\phi[y]f_2[x-y]\\
&=aH\{f_1[x]\}+bH\{f_2[x]\}\end{align}$$

So, the linearity holds. There is nothing special about the fact that the system is 2-dimensional or similar. You just follow the rules for the calculation, and you end up with the statement of linearity.

A similar thing can be carried out for the shift-invariance.


Computational complexity lies in the correlator banks you need to despread the GPS signal, not in the position and time calculation (these are not separate things. You can't have one without the other). In fact, the TDOA system is very likely far more computationally intense, since it needs to do more cross-correlations.

You need to use the first method; 2. can't work, you'd need a clock that is about as stable and thus as large, power-hungry and expensive as that on the GPS satellites. Even with that clock, it wouldn't work well, because you still need to account for a lot of effects, including special relativity, that you can't solve via simple TDOA.

Computational complexity is still rather high, indeed, but really, not a problem since the late 1990s, where highly integrated chips for GPS became the norm. What it does do is drain the battery of your smartphone!


QPSK modulation is a PSK modulation with exactly 4 constellation points.  The constellation points for QPSK are normally at $\left(\pm\dfrac{1}{\sqrt{2}}\pm\dfrac{1}{\sqrt{2}}j\right) $


How can I shape my time channel response, in order to have Gaussian shaped Doppler Spread?

See  section 2.4.3 of this reference https://web.stanford.edu/~dntse/Chapters_PDF/Fundamentals_Wireless_Communication_chapter2.pdf

If the doppler spectrum has to be gaussian, the auto-correlation of tap gains should be gaussian (which is correctly mentioned in other answer but I somehow felt more details were not captured). For a coherence time $T_c$, proportional to doppler spread $1/D_s$, will be expressed as smallest amount of $n$ for which your tap at $l$ will be significantly different from tap at $l+n$. So for example, if $n=5$, gain of taps from $0$ till $4$ will be having a guassian distribution, while the gain at tap $0$ will be uncorrelated with tap at $5$.

If you model your taps as Rayleigh fading channel, your PSD will not be Gaussian as given in eq 2.60 of the reference. You may need to generate your tap gains such that their auto-correlation function follows gaussian distribution.

Since you are using MATLAB

raychan = comm.RayleighChannel('DopplerSpectrum',{doppler('Gaussian')});

This makes sure that each of the channel taps gain auto-correlation (and PSD) follows Gaussian distribution.


That's one way to do it but but typically you would have a single line of three delays in the feedback path and apply $a$ after the first delay, $b$ after the second and $c$ after the third.

And yes, you need arrowheads to make this clear.

By convention the coefficient of the feedback path would be called $a_1$, $a_2$, $a_3$ and the feed forward path would be $b_0$, $b_1$, etc.

See for example https://www.researchgate.net/figure/Direct-Form-I-structure-of-IIR-filter_fig1_236688224


Gaussian wavelet generation of a given frequency

The modulation property of the Fourier Transform states that given a continuous time domain function $x(t)$ with Fourier transform $X(f)$, the function
$h(t) = e^{j 2\pi f_0 t} x(t)$ has the Fourier transform $X(f-f_0)$ (i.e. shifted by $f_0$ in the frequency domain):

$$
\begin{align}
   e^{j 2\pi f_0 t} x(t) \overset{\mathcal F}\Longleftrightarrow X(f-f_0)
\end{align}
$$

Similarly for a discrete time signal $x[n]$ with period $N$ and Discrete Fourier Transform $X_k$ (under an $N$-point DFT), the equivalent relation is:

$$
\begin{align}
   e^{j 2\pi \frac{m}{N} n} x[n] \overset{\mbox{DFT}_N}\Longleftrightarrow X_{k-m}
\end{align}
$$

As you might have noticed, these expressions multiply you original time-domain signal by a complex term, which gives you a complex valued time-domain signal. If you want a real valued time-domain signal, you would actually need to construct a combination of a signal frequency shifted by $m$ with another one frequency shifted by $-m$. This can be carried out in the time-domain by multiplying the signal by $\cos(2\pi\frac{m}{N} n)$.


Why are males' pitch lower than females?

Males usually have deeper voices than females because they have longer vocal folds.

Yes

Short and relaxed vocal folds produce a deep voice, long and tense vocal folds produce a high-pitched voice.

Yes on "relaxed", no on "short".
Treating the vocal cords as strings: The resonance frequency is a function of mass (or density), tension and length.
Higher mass -> lower pitch
Higher tension -> higher pitch
Higher length -> lower pitch

Male vocal cords are longer and heavier.


Finding when a part of audio starts and ends in another audio

You may consider using Dynamic Time Warping (DTW) for this.
One approach which generally would work is as follows:

Divide both the audios into smaller (analysis) windows with appropriate overlap as per the needed temporal resolution
Transform the audio from time domain to cepstral. You may compute  MFCCs.
Copute DTW to analyse the warping


For the number of samples given, your initial statement is false.  Sampling at exactly T/2 is a mathematical limit.

For any finite number of samples, you need to sample a band-limited pure sinusoid of period T more frequently than T/2 to avoid aliasing or other pathologies, or at greater than 2X the highest spectrum frequency (for a baseband range). (and in real-life, to take into account the transition bandwidth and stop-band of any low-pass anti-aliasing filtering).  The shorter the sampling interval, the higher the sample rate needs to be above 2*Fmax (for a given baseband band-limited range) to reduce numerical, quantization, and other signal-to-noise induced errors to below some given threshold.

As you approach an infinite number of samples (beyond the life of the universe), you approach being able to reconstruct a sinewave with sampling at T/2.


What is the physical significance of negative frequencies?

Negative frequency doesn't make much sense for sinusoids, but the Fourier transform doesn't break up a signal into sinusoids, it breaks it up into complex exponentials (also called "complex sinusoids" or "cisoids"):

$$F(\omega) = \int_{-\infty}^{\infty} f(t) \color{Red}{e^{- j\omega t}}\,dt$$

These are actually spirals, spinning around in the complex plane:

(Source: Richard Lyons)

Spirals can be either left-handed or right-handed (rotating clockwise or counterclockwise), which is where the concept of negative frequency comes from.  You can also think of it as the phase angle going forward or backward in time.

In the case of real signals, there are always two equal-amplitude complex exponentials, rotating in opposite directions, so that their real parts combine and imaginary parts cancel out, leaving only a real sinusoid as the result.  This is why the spectrum of a sine wave always has 2 spikes, one positive frequency and one negative.  Depending on the phase of the two spirals, they could cancel out, leaving a purely real sine wave, or a real cosine wave, or a purely imaginary sine wave, etc.

The negative and positive frequency components are both necessary to produce the real signal, but if you already know that it's a real signal, the other side of the spectrum doesn't provide any extra information, so it's often hand-waved and ignored.  For the general case of complex signals, you need to know both sides of the frequency spectrum.


The "wavespeed" is the speed of light, which is 299,792,458 m/s if in free-space, or slowed by the root of the relative permittivity (dielectric constant for the material where the wave is propagating).  So $v = c/\sqrt{\epsilon_r}$ and the formula for wavelength is $\lambda = c/(f\sqrt{\epsilon_r})$, where $c$ is the speed of light (299,792,458 m/s), $f$ is the frequency of the signal in Hz, and $\epsilon_r$ is the dielectric constant.

Example dielectric constants:

Air at 0°C: 1.000594
Air at 20°C: 1.000528
Glass: 5 to 10
Paper: 3.85
Silicon: 11.68

For example, a 30 GHz signal in air would have a wavelength of nearly 10 mm: $ 299,792,458/30\times10^{9}$.