What level of correlation coefficient suggests correlation?

One way to look at it is that $r^2$ (yes, I do mean $r^2$, not $r$ or $|r|$) is the fraction of the variance of the target (the series A) that is explained by the estimator (a linear function $\alpha + \beta B$ of the series B).  Note that if Series A had variance $\sigma^2$, then the variance of $(A - (\alpha + \beta B))^2$ (what B cannot explain about A) is $(1-r^2)\sigma^2$.  So, if $|r| < \frac{1}{\sqrt{2}} \approx 0.707\ldots$, less than half the variance has been explained, and so one viewpoint is that a correlation of less than $70\%$ or so suggests that the linear estimator being considered here is not all that great.


What kind of filter would have a sawtooth frequency response?

A windowed design is easy. Starting with an ideal zero-phase frequency response:

$$\frac{\omega}{2\pi} + \frac{1}{2}$$

The inverse Fourier transform of that gives the ideal impulse response:

$$\int^\pi_{-\pi}\left(\frac{\omega}{2\pi}+\frac{1}{2}\right)e^{i\omega k}d\omega = \frac{\sin(\pi k)}{k} + i\left(\frac{\sin(\pi k)}{\pi k^2} - \frac{\cos(\pi k)}{k}\right).$$

Similar to the $\text{sinc}$ function, the real part $\frac{\sin(\pi k)}{k}$ has an illegal division by zero at $k=0$ where the limiting value $\pi$ should be used instead. For other integer $k$, the real part is zero. The imaginary part has similar problems at $k=0$, where its limit is 0. Drawing from Rick's answer, at integer $k$, $\sin(\pi k)$ is zero and $\cos(\pi k)$ is +1 or -1. This simplifies the ideal impulse response at integer $k$ to:

$$\begin{cases}\pi&\text{if }k=0,
\\i\frac{-(-1)^k}{k}&\text{otherwise.}\end{cases}$$

If the filter is implemented as a sum of a purely real and a purely imaginary filter, the purely real filter requires no work other than to multiply the input by a constant. The imaginary part of the impulse response is antisymmetrical, which may enable savings.

To continue with the windowed design, choose a window function and its length (equals the number of filter taps) and multiply the ideal impulse response by it.
That gives the impulse response of the filter. The window function should be centered at $k = 0$.


Does FFT magnitude of recorded sound depends on the bit depth used?

With very few exceptions any audio processing these days is done in floating point, i.e. the data that comes through the ADC is immediately converted to floating point numbers and normalized so that the max voltage corresponds to a value of 1.0. This eliminates any dependency of the scaling to the number of bits.

Fixed point processing is quite difficult and typically requires specialized data types and intensive scaling and noise management.

The output of an 8-bit ADC could be interpreted as the numbers $-128,-127, ...+127$ but it's almost always better to think of it as $-1, -127/128, ..  -1/128 , 0 ,+1/128 ... +127/128$


Understanding the definition of signal's bandwidth

"tutor explained that the bandwidth indicates how fast the signal is changing"

This is not correct even in a hand-waving fashion.  As you were told, a pure sinusoid $x(t) = \sin(2\pi f_0t + \theta)$ has zero bandwidth (according to the definition of bandwidth that your tutor is using) but that does not mean that the sinusoidal is not changing at all! It is changing and its maximum slope (a.k.a. "how fast the signal changing") is
$$\max \left| \frac{\mathrm dx}{\mathrm dt}\right| = \max \left|2\pi f_0 \cos(2\pi f_0t+\theta)\right| = 2\pi f_0.$$ More generally, a result due to Bernstein says that if the maximum frequency in a generic $x(t)$ bounded within $[-1,1]$ is $f_0$, that is, $X(f) = 0$ for $|f| > f_0$, then
$$\max \left| \frac{\mathrm dx}{\mathrm dt}\right| \leq 2\pi f_0.$$
In short, the minimum frequency (in a bandpass signal) is irrelevant; it is the maximum frequency that matters, and the "how fast the signal is changing" depends on the maximum frequency and not on the bandwidth.

Turning to the rectangular pulse that puzzles you, note that the signal doesn't have a derivative (as far as mathematicians, including Bernstein, are concerned) when the signal makes the abrupt transition from high to low or vice versa but here on dsp.SE we have no qualms about saying that the derivative is $\pm \infty$ at these points.  As MattL says in a comment, if you know the Fourier transform of a rectangular pulse (hint: $\text{rect}(t) \leftrightarrow \text{sinc}(f)$), then you know that the rectangular pulse has no maximum frequency (there is no finite $f_0$ such that $\text{sinc}(f) = 0$ for all $|f|>f_0$) and that the infinitely fast changing of the rectangular pulse is well-supported by the infinite bandwidth of the rectangular pulse.


Finding the friction coefficient by using experimental data

I can only propose a crude starting point based on a more modest work I did. My goal was to determine the bleaching coefficient of a fluorophore when submitted to varying illumination intensity over time. I had the theoretical law for the bleaching over time as a function of intensity (depending on unknown parameters), the amount of fluorophore over time, and the starting intensity and intensity at the end of the experiment. The goal was to retrieve the parameters for the bleaching law and the whole intensity over time.

Anyway. What I did was to construct (using MATLAB) a function to minimize. In your case, given that you know $x(t)$ experimentally, I suggest to build a function $O(x, y)$ such that:

$$ O(x, y) = \sum_t |x(t) - y(t)|^2 $$

with $x(t)$ being your experimental data and $y(t)$ being the solution of your differential system for some parameters $k, M, a$. Put this function as the objective in a simplex optimizer that will find its minimum, allowing $k, M, a$ to vary.

This is done through numerical solution. So you have a function to minimize, say $F(k,a,M)$ that returns the value of $O(x,y)$ (which is scalar btw) for the current numerical resolution with the given $k,a,M$.

At each step, the whole system will be solved again, and the solution will be matched against the experimental data until they fit.

I have found this to be working quite well in my case, as long as you are sure you have the correct model to describe the experimental data. Otherwise, expect some seriously irrelevant results.


Prerequisites for audio signal processing

DSP is indeed a very mathematical discipline, but not completely. The amount of mathematical knowledge you need to operate DSP concepts comfortably is luckily limited to a particular fairly small subset. I would say that you need the following:

Complex variables. DSP deals with oscillating signals and systems, and those are very conveniently represented using Euler's complex exponential formulas. This includes knowledge of polynomials over complex variables, and so on.
Simple differential calculus. I would say that you need to know how to take derivatives of polynomials and other common functions such as exponentials and logarithms. Knowing l'Hôpital's rule also occasionally comes in handy. Taking derivatives is primarily required for finding minima and maxima of functions. Note that you need to be able to do this for complex variables as well, but at this stage there is practically no difference in concepts between real and complex differentiation.
Linear algebra. I put linear algebra before integration because it is much more useful, and since integration is a linear operation, a lot of it can be replaced by linear algebra concepts. I can't exaggerate how important linear algebra is, so go and learn it. I highly recommend Gilbert's MIT Open Courseware linear algebra class. No one teaches it better than this guy. Just watch the lectures.
Integral calculus. I'm not quite sure how to approach describing the importance of integration in DSP. If linear algebra was a moderately easy, fun and useful topic to learn, integration, especially complex integration, takes a lot of commitment and intuition, both of which require lots of time. Even though I usually advocate for in-depth mathematical knowledge, here's what you can very easily get away with. Learn about integration concepts, what it means and how integration is related to derivatives. Since most integrals are difficult to work out by hand, people use integration tables to look up answers instead. This will get you through most integration you need to know for DSP.
Infinite series. This is the one part of calculus that everyone forgets the minute they take that final exam, but it actually comes in handy in DSP. Specifically, when analyzing recursive systems, i.e. those that take their output and reuse it as input creating a feedback loop, you need to be able to spot several types of infinite number patterns that appear and point to a single formula that produces every single term in that infinite string of numbers or polynomial coefficients.
Probability. I hesitate a little to write about this one, because not all DSP courses strictly require probability, but if the one you're looking at does, it wouldn't hurt to know it. Not be be lost, you need to know expressions and concepts of means and variances, as well as Bayes' rule.

Obviously, what I've given above are the minimal requirements, but you may need to learn more as you go.


Given Position Measurements, How to Estimate Velocity and Acceleration

One approach would be to cast the problem as least-squares smoothing. The idea is to locally fit a polynomial with a moving window, then evaluate the derivative of the polynomial. This answer about Savitzky-Golay filtering has some theoretical background on how it works for nonuniform sampling.

In this case, code is probably more illuminating as to the benefits/limitations of the technique. The following numpy script will calculate the velocity and acceleration of a given position signal based on two parameters: 1) the size of the smoothing window, and 2) the order of the local polynomial approximation.

# Example Usage:
# python sg.py position.dat 7 2

import math
import sys

import numpy as np
import numpy.linalg
import pylab as py

def sg_filter(x, m, k=0):
    """
    x = Vector of sample times
    m = Order of the smoothing polynomial
    k = Which derivative
    """
    mid = len(x) / 2
    a = x - x[mid]
    expa = lambda x: map(lambda i: i**x, a)
    A = np.r_[map(expa, range(0,m+1))].transpose()
    Ai = np.linalg.pinv(A)

    return Ai[k]

def smooth(x, y, size=5, order=2, deriv=0):

    if deriv > order:
        raise Exception, "deriv must be <= order"

    n = len(x)
    m = size

    result = np.zeros(n)

    for i in xrange(m, n-m):
        start, end = i - m, i + m + 1
        f = sg_filter(x[start:end], order, deriv)
        result[i] = np.dot(f, y[start:end])

    if deriv > 1:
        result *= math.factorial(deriv)

    return result

def plot(t, plots):
    n = len(plots)

    for i in range(0,n):
        label, data = plots[i]

        plt = py.subplot(n, 1, i+1)
        plt.tick_params(labelsize=8)
        py.grid()
        py.xlim([t[0], t[-1]])
        py.ylabel(label)

        py.plot(t, data, 'k-')

    py.xlabel("Time")

def create_figure(size, order):
    fig = py.figure(figsize=(8,6))
    nth = 'th'
    if order < 4:
        nth = ['st','nd','rd','th'][order-1]

    title = "%s point smoothing" % size
    title += ", %d%s degree polynomial" % (order, nth)

    fig.text(.5, .92, title,
             horizontalalignment='center')

def load(name):
    f = open(name)
    dat = [map(float, x.split(' ')) for x in f]
    f.close()

    xs = [x[0] for x in dat]
    ys = [x[1] for x in dat]

    return np.array(xs), np.array(ys)

def plot_results(data, size, order):
    t, pos = load(data)
    params = (t, pos, size, order)

    plots = [
        ["Position",     pos],
        ["Velocity",     smooth(*params, deriv=1)],
        ["Acceleration", smooth(*params, deriv=2)]
    ]

    create_figure(size, order)
    plot(t, plots)

if __name__ == '__main__':
    data = sys.argv[1]
    size = int(sys.argv[2])
    order = int(sys.argv[3])

    plot_results(data, size, order)
    py.show()

Here are some example plots (using the data you provided) for various parameters.

Note how the piecewise constant nature of the acceleration becomes less obvious as the window size increases, but can be recovered to some extent by using higher-order polynomials. Of course, other options involve applying a first derivative filter twice (possibly of different orders). Another thing that should be obvious is how this type of Savitzky-Golay filtering, since it uses the midpoint of the window, truncates the ends of the smoothed data more and more as the window size increases. There are various ways to solve that problem, but one of the better ones is described in the following paper:

  P.A. Gorry, General least-squares smoothing and differentiation by the
  convolution (Savitzky–Golay) method, Anal. Chem. 62 (1990) 570–573. (google)

Another paper by the same author describes a more efficient way to smooth nonuniform data than the straightforward method in the example code:

  P.A. Gorry, General least-squares smoothing and differentiation of
  nonuniformly spaced data by the convolution method, Anal. Chem. 63
  (1991) 534–536. (google)

Finally, one more paper worth reading in this area is by Persson and Strang:

  P.O. Persson, G. Strang, Smoothing by Savitzky–Golay and Legendre
  Filters, Comm. Comp. Finance 13 (2003) 301–316. (pdf link)

It contains much more background theory, and concentrates on error analysis for choosing a window size.


How do audiobook players manage to quickly speedup/slowdown playback speed without using more memory?

WOW long gap to answer this question, maybe if I have a Einstein brain I can understand the image, what is not my case lol.
I haved touched in this algorithm in a long time ago, and I remember that was a surprised for me that almost no one around the Web talks about the Algorithm from Morita Naotaka, I tryed seach his original paper but I cant found online, I think that it is in Japanese, well I dont know ...
So I learned about it looking at among number of patents that make references to the Morita Naotaka, and I freeze when I see a lot of patents around 2000 years that use the same code.
As far as I know this algorithm was derived from TDHS, has advantages about TDHS, this does the job using just Period size Buffer. And believe me it's as easy to code as TDHS.
So lets go, take a look at the Yamaha Patent number US8085953 this is basicly the PICOLA from Morita Naotaka...
Here one image from this Patent that is infinity better to understand how this works to compress :

Looks like one TDHS using a triangular window, but take a look how it are cutting wave A and wave B in the exact Period size(Lp), then they are crossfaded A + B. this image seems to be an compression using rate = 0.6 factor(alpha),  L can be determined using:
$$ \begin{align}
L & = \operatorname{round}\left(\frac{L_p}{\alpha-1}\right)
\end{align} $$
And here how its works to expansion:

Now L can be determined using:
$$ \begin{align}
L & = \operatorname{round}\left(L_p\frac{\alpha}{1-\alpha}\right)
\end{align} $$
Time ago I show a basic code to compress based in TDHS here, now lets try a litle bit more complex and complete example to expand (it will slowdown), so based in this images/equations, here my PICOLA Algorithm:
alpha=0.5;

f=735;
Fs=44100;
signal= 0.9*sin(2*pi*f/Fs*(1:44100)); signal = signal';

period= 60; % is it (Fs/F)

[signal,Fs]=audioread('ederwander.wav');

nsamples = length(signal);
out=zeros(1,floor(nsamples/alpha));
inptr=1;
outptr=1;

minF = 50;
maxF = 900;

minP = floor(Fs/maxF);
maxP = floor(Fs/minF);

winsize = 1024;

while ( inptr+winsize <= nsamples )

    %chunk used to find the period, its need be >= MaxP
    windowedSignal=signal(inptr:inptr+winsize);

    %Basic Autocorrelation to find the period
    Maximum=-Inf(1);
    for P =minP:maxP
        ac = sum(windowedSignal(1:(winsize)-P) .* windowedSignal(P:(winsize -1)));
        if ac > Maximum
            period=P;
            Maximum=ac;
        end
    end

    period = period-1;

    %split all how shown in the patent picture
    waveA = signal(inptr:inptr+period-1);
    waveB = signal(inptr+period:(inptr+period*2-1));
    CrossfadeAB = (waveA .* linspace(0,1,period)') + (waveB .* (1-linspace(0,1,period))');
    L = round(period*alpha/(1-alpha));
    waveB_T0 = signal(inptr+period:inptr+L);
    out(outptr:outptr+period-1)=waveA;
    out(outptr+period:outptr+period*2-1)=CrossfadeAB';
    out(outptr+period*2:outptr+period+L)=waveB_T0';

    %increment input and output pointers
    inptr = inptr + L;
    outptr = outptr + L + period;

end

%not best choice, just append the end to match the out size, this will click, maybe need crossfade here
if outptr < nsamples/alpha

    fim=floor(nsamples/alpha)-outptr;
    out(outptr:floor(nsamples/alpha))=signal(nsamples-fim:nsamples);

end

plot(out);
sound(out,Fs);

I think that standard algorithm of PICOLA just get nice results to 0.5<=alpha<=2
And one more time it's can be considered a proof of concept to show how its can be done!


Convert normalized frequency to real frequency in AR model

Matlab generally uses normalized frequencies in radians/sample. If the signal is real-valued then these frequencies are in the interval $[0,\pi]$ (as you can see in the plot). The Matlab frequency variable $w$ is normalized as

$$w=\frac{2\pi f}{f_s}$$

where $f_s$ is the sampling frequency. So $w=\pi$ corresponds to $f=f_s/2$.


What Are the Differences between Super Resolution, Denoising and Deblurring?

All 3 of them fall into the category of Inverse Problems in the Image Processing world.

Lets assume a Linear Model and then we will be able to show all 3 of them as parameters of the same framework. Then the differences will be clear and one could generalize it into Non Linear settings as well.

Have a look on the following model:

$$ y \approx A x + w $$

Where $ x $ is the image we want to infer from the given image $ y $, $ A $ is the degradation operator applied on $ x $ and $ w $ is white additive noise.

Denoising
For denoising only noise is added hence $ A = I $ and we left with estimating $ x $ from a noisy measurements.
Deblurring
In that case $ A $ is a matrix form of some Low Pass Filter (Circulant Square Matrix) which applies a blur on the image. We try to infer $ x $ from a blurred and noisy version of it given by $ y $.
Super Resolution
In this case take $ A $ from (2) and remove few rows of it. Namely we get only some of the data in $ x $ and the data is both blurred and noisy and we try to infer the whole data of $ x $ given $ y $.

Usually the problems above aren't well defined which means they require some kind of regularization.

Have a look on the following problem:

$$ \arg \min_{x} \frac{1}{2} {\left\| A x - y \right\|}_{2}^{2} + \lambda r \left( x \right) $$

This is the Linear Inverse problem where you can adapt to any of the cases above. The $ r \left( \cdot \right) $ is the regularization function which means we have some model on the data. For instance for Images many use the Total Variation model which means the Image is Piece Wise Smooth (Sparse Derivatives).


Difference between Fourier Transform and DFT? - Example

The difference is between $\mbox{sinc}$ and the periodic version obtained using the DFT.

See this answer for a comparison.

  It strikes me that asinc = ratio of two sincs.

The $\mbox{asinc}$ function is a ratio of two $\sin$ functions:
$$
\frac{\sin(N\omega/2)}{\sin(\omega/2)} = \frac{\sin(N\omega/2)}{\omega/2} \frac{\omega/2}{\sin(\omega/2)}
$$
so, yes, you could see it as the ratio of two $\mbox{sinc}$ functions.

  Underlying question: can I model the sampling of a continuous signal (by, say, a digital camera) in the frequency domain as the DFT of the signal convolved with a Dirac delta train times the DFT of the pixel footprint (a rect)?  If so, why is the DFT of the pixel footprint an asinc instead of a sinc?

So you seem to be trying to model taking a picture with a digital camera?   Camera modeling is a little tricky. Even Gonzales and Woods seem to gloss over it.

The fundamental issue is the DFT of a $\mbox{rect}$ ($\Pi$) is a $\mbox{asinc}$.

If you're doing a discrete-time Fourier transform (DTFT), then it's not, but usually when dealing with computed FTs, you want the DFT.


Main thing is that in the first frame it will be required to select the object of detection because it is obvious that the algorithm will not know automatically which object you want to track if a number of things will be moving in the video scene.Lucas-Kannade method is one of the methods which can detect moving objects in a given video frame . If you do not specify object it will show movement of all moving objects.Thats why you have to tell the algorithm which thing you want to track. Second thing i want to tell you is about Kalman filter which i think you are interpreting wrong , what Kalman filter actually do is that it uses the system model and input signal and generates its own output which is actually a predicted value . Now using the observation from the system it calculate the correction to be applied to the predicted output and hence final results are better than observations obtained from the system. Hence very first step is that you need a system model and the measurements. In tracking algorithms generally, Image processing algorithm is applied to detect the object of interest and then Kalman filter is applied to make the detected location more precise. Do not think that it will predict on its own. I hope it will help you.


With typical sampling (not noise shaped) you get half a bit for every halving of the sampling rate (with proper filtering of the out of band noise prior to downsampling). So a decimate by two for example consists of a half band filter followed by a down-sample by two which is done simply by selecting every other sample.  Assuming white noise, the half band filter removed half the noise or 3 dB, and the the SNR of ADC’s goes as 6 dB/bit.
Therefore, if you went from 30 MHz to 7.5 MHz with a decimate by four, you would gain a full bit or 12 bits extended to 13 bits assuming noise is limited by quantization noise and the quantization noise process is sufficiently white as typically assumed and approximated.
It’s should be straight forward to see now how to go from 12 bits to 14 bits through successive decimations. (And the filtering and decimation can be done in one block each or distributed which gives flexibility in the filter designs used).
An additional caution when working with real hardware is to refer to the expected Effective Number of Bits (ENOB) for the actual ADC used and in what conditions of sampling, input frequency etc as the number of bits to be achieved and improved on, and pay attention to actual spur levels (as given by the Spurious Free Dynamic Range, SFDR) that will ultimately limit the SNR since decimation will do nothing to reduce the power level of in band spurs. Often these spurs are below the noise level due to quantization in which case dynamic range improvement can be gained through decimation.


Sensor fusion with accelerometer and gyroscope only

Actually, the plots for the components of gravity look good. What you miss is the fact, that your overall gravity (i.e. gravity sum) is the length of the vector, that is spanned by the three components. Hence, you should calculate

$$G=\sqrt{g_x^2+g_y^2+g_z^2},$$

where $g_x,g_y,g_z$ are the three single components measured by the acelerometer. Looking at your curves, I predict the sum will be quite flat.


Manually Windowing an IIR Filter in Matlab

In the classical world of digital filters, an infinite impulse response (IIR) filter is used for its computational efficiency while achieving supreme performance compared to a finite impulse response (FIR) filter.

Usually, implementing an IIR filter happens in the time domain, with an architecture of direct form-I, II or else which provides a solution of the linear constant coefficient difference equation (LCCDE) describing the filter.

$$
\sum_{k=0}^{N} a[k]y[n-k] = \sum_{k=0}^{M} b[k]x[n-k]
$$

which is typically called an $N^{th}$ order LCCDE.

However, when there happens a need for implementing the IIR filter using time domain convolution instead;

$$
y[n] = h[n] \star x[n]
$$

then this requires the practical availability of the impulse response $h[n]$ of the IIR filter which is infinite length by definition.

To obtain a finite length, FIR, approximation $h[n]$ to the IIR digital filter, first choose your analog prototype filter, then apply the necessary transforms to get the digital IIR equivalent filter coefficients $a[k]$ and $b[k]$, and compute the truncated impulse response $h[n]$ of this filter, from $a[k]$ and $b[k]$ to a length that's adequate for the purpose of application and finally apply a further windowing to $h[n]$ if that's necessary...

The following MATLAB/OCTAVE command will give you a length $N$ truncated impulse response from IIR coefficients $a[k]$ and $b[k]$ with zero initial conditions.

h = filter(b,a,[1 zeros(1,N-1)]);


differential phase tracking loop for the receiver

There is no need to track the phase difference at all.  What is needed is to compare the "average" RF phase during the current bit interval to the "average" RF phase during the previous bit interval.  If the two phases are the same, or near enough the same, the receiver decides that the phase did not change during the transition from the previous bit interval to the current bit interval and so the differential encoding at the transmitter is telling us that a $0$ data bit is being transmitted.  If the two phases are vastly different, then we are being told that a $1$ data bit is being transmitted.

How is all this done? Well, use the part of a QPSK receiver that separates the RF signal into I and Q signals, and then match-filter (or otherwise process) the I and Q signals to produce complex-valued outputs $Z_1 = (X_1,Y_1)$ and $Z_2 =(X_2,Y_2)$
where $X$ is
the output of the I branch and $Y$ the output of the Q branch of the receiver. Note that the $X$'s and $Y$'s are individual real numbers: they are the sample values obtained from the I and Q branch filters/processors at the chosen sampling instants, not continuous-time signals or vectors of sample values that DSP processors or computers use to represent continuous-time signals.

The  hard-decision DBPSK decision device considers the question:

  Is the new symbol $Z_2$ closer to the old symbol $Z_1$ or closer to the
  negative $-Z_1$ of the old symbol?

and thus compares

$$(X_2-X_1)^2 + (Y_2-Y_1)^2 \gtrless (X_2+X_1)^2 + (Y_2+Y_1)^2$$

which can be simplified to a sign comparison on $\langle Z_1,Z_2\rangle = X_1X_2+Y_1Y_2$.
Note that this is essentially asking

  Are the two vectors $Z_1$ and $Z_2$
  are pointing in roughly the same direction (in which case the inner product
  or dot product is positive) or in roughly opposite direction (in which case
  the dot product is negative)?

A third viewpoint thinks of $Z_1$ and $Z_2$ as
complex numbers and asks

  Is $\text{Re}(Z_1Z_2^*) = X_1X_2+Y_1Y_2$ positive or negative?

But, regardless of the viewpoint, the end result is the same.

Note that there is no explicit phase comparison in all this: those boxes marked $\arctan\left(\frac yx\right)$ on block diagrams of communications receivers are more for conceptual understanding than for implementation guidance.

Some notes to read before downvoting this answer.

It is assumed that the receiver has a good estimate of the received carrier frequency but not of the carrier phase.  Frequency acquisition can be done by squaring the RF signal (to eliminate the phase modulation), high-pass filtering (to remove the DC in the squared signal which has DC plus an unmodulated double-frequency signal in it), and frequency-halving, to recover the carrier frequency. In fact, the output of the frequency halver can be used as the local oscillator in the receiver front end.  It is not necessary to track the RF phase.  All that we are interested in is whether the phase in the current interval was the same as the phase in the previous interval or was it different? Where exactly $Z_1$ and $Z_2$ individually lie in the complex plane is not of importance; their relative position is.  This is different from phase-coherent demodulation of plain vanilla BPSK in which there is a single (real-number) output that must be compared to $0$ to determine the transmitted bit and woe betide the designer who ignores the possibility of phase ambiguity caused by the frequency halver.  Note also that we are assuming that the bit interval timing is known with a fair degree of accuracy so that the DBPSK receiver knows when it is the correct time to sample the I and Q branch outputs to get the $Z$'s.


Closed form expression for $\sum_{k=0}^{n} \alpha^{-k}$?


Why do poles in the left half of the S plane make a system stable?

An $n$th order linear system is asymptotically stable only if all of the components in the homogeneous response from a finite set of initial conditions decay to zero as time increases, or mathematically written:
$\displaystyle{\lim_{t \to \infty}}\sum_{i=1}^{n}C_ie^{p_it}=0$ ($p_i$ are the system poles).  As time increases (in a stable system) all components of the homogeneous response must decay to zero.If any pole has a positive real part there is a component in the output that increases without bound, causing the system to be unstable. So, in order for a linear system to be stable, all of its poles must have negative real parts (they must all lie within the left-half of the s-plane).  An "unstable" pole, lying in the right half of the s-plane, generates a component in the system homogeneous response that increases without bound from any finite initial conditions.


Inverse $\mathcal Z$-transform when region of convergence goes outwards from the inner pole?

There exist $3$ sequences with the given expression as their $\mathcal{Z}$-transform. These $3$ sequences correspond to $3$ different regions of convergence (ROCs):

$|z|>\frac12$: right-sided
$\frac13<|z|<\frac12$: two-sided
$|z|<\frac13$: left-sided

It's important that you figure out both right-sided and left-sided inverse $\mathcal{Z}$-transforms of the basic transform

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

and their respective ROCs. Then it's easy to find all $3$ inverse transforms for the given example.


modelling wireless multi-path channel - is channel response sum of delta(t-path_delay) or "delta" should be smoothed somehow?

Depends on the physics of the situation. The model is reasonably good if you have secular reflections and it is not very good if you have diffuse reflections. If the size and geometry details of the reflectors are large as compared to the wavelength, than its mostly specular. If the surfaces are uneven are very non flat, then its more diffuse. This is all measured in units of wavelengths. For a 5GHz link that's about 6cm. For an FM radio link that's about 3m.
It also depends on the relative bandwidth of the signal. For most wireless links that's pretty narrow, i.e the behavior at the carrier frequency is a good model for the behavior of the overall system. Most natural physical transfer functions such as a reflection of a wall do not change very rapidly with frequency so that's why the singe impulse is a good starting point.
Of course there are always exceptions to the rules so it depends on your specific problem.


Estimate the wavelength of an oscillatory signal

Consider the code snippet below to find the three peaks in the graph. Note, that peak search usually only searches for positive peaks. Hence, I negated the signal to search for the down-peaks. Also, I adapt the SNR threshold in the code to detect only the very strong peaks.

import urllib2
import scipy.signal as signal

url = 'https://raw.githubusercontent.com/Omer80/wavelength/master/oscillations.dat'
response = urllib2.urlopen(url)
cr = np.loadtxt(response)
X = cr[0,:]
Y = cr[1,:]

plt.subplot(2,2,1)
plt.plot(Y)

plt.subplot(2,2,2)
plt.plot(X, Y)

peakind = signal.find_peaks_cwt(-Y, np.array([100]), min_snr=5)
peakind = np.array(peakind)
print peakind

plt.stem(X[peakind], 1+0*peakind);

period = int(np.mean(np.diff(peakind)))
print "Period between peaks:", X[period]-X[0]

Output:

[1919 2500 3080]
Period between peaks: 0.145036259065


sketching frequency domain representation

This is a classical example of beat. This means that the signal (assuming that it is periodic and extends to plus and minus infinity) contains only two frequencies with the same amplitude. The frequency of the sine wave which seems to enclose the signal will be equal to the difference between the two frequencies. The frequency of the sine wave which seems to be modulated by the enclosing wave will be the mean of the two frequencies which are present in the signal.


Find Bandwith from Spectral Efficiency and Eb/N0

Your initial answer is correct -- the bandwidth depends only on the spectral efficiency and the symbol rate. $E_b/N_0$ is a figure of merit largely designed to be independent of the bandwidth.

Let me give you an example of where $E_b/N_0$ might come into play. Say you are required to provide $E_b/N_0 = 10$, and you know that your receiver is rated at $N_0 = 10 ^ {-12}$. Then, $E_b = 10 N_0 = 10^{-11}$. At a rate $R = 2 \times 10^6$, the received power is $(2 \times 10^6)(10^{-11} = 2 \times 10^{-5}$ watts, or $-47$ dB.

If you know the channel attenuation, then you can calculate the transmitted power. From there, you can calculate things like battery life. And so on...


Concept of combining multiple FIR Filters into 1 FIR filter

What you tried to do in your example is combine the filters in series rather than parallel.
That is, if your signal is $$y = ((x*h_1)*h_2)\ldots*h_7$$ then, according to the associativity property of the Convolution, $$y=x*\underset{H}{\underbrace{(h_1*h_2\ldots*h_7)}}$$

However, if you wish to perform a parallel filtering then your signal is $$y=\sum_{i=1}^7 x*h_i$$
and acording to the distributivity property of the Convolution, $$y=x*\underset{H}{\underbrace{\sum_{i=1}^7 h_i}}$$


So I think it's mostly missing the line

    dft2d = np.zeros((M,N),dtype=complex)

in your code. Otherwise your sum may not end up being complex.

And then you can't then plot the resulting image... because that needs to be real-valued.

I modified your code so that it sort of works (but is very slow). See below.  This results in something more like what you expect --- though because of the indexing, the peak is in the top-left corner rather than the middle.  There's probably something like Matlab's fftshift function that will center it for you.

Updated Code Below

import numpy as np
import matplotlib.pyplot as plt
import PIL
import cmath

def DFT2D(image):
    data = np.asarray(image)
    M, N = image.size # (img x, img y)
    dft2d = np.zeros((M,N),dtype=complex)
    for k in range(M):
        for l in range(N):
            sum_matrix = 0.0
            for m in range(M):
                for n in range(N):
                    e = cmath.exp(- 2j * np.pi * ((k * m) / M + (l * n) / N))
                    sum_matrix +=  data[m,n,1] * e
            dft2d[k,l] = sum_matrix
    return dft2d

img = PIL.Image.open("example.jpg")
img2 = img.resize((50,50))
plt.imshow(img2)
dft = DFT2D(img2)
plt.imshow(dft.real)


Using the Euler approximation, the mapping between the $s$-plane and the $z$-plane would be:

$$
  s \rightarrow \frac{z - 1}{T}
$$
where $1/T = f_s = 30$ (thus $s \rightarrow 30(z-1)$).

Substituting back in the original transfer function yields:

$$
\begin{align}
  H(z) &= \frac{24.2\cdot 30 (z-1) + 1}{2.2\cdot 30 (z-1) + 1} \\
       &= \frac{726 (z-1) + 1}{66 (z-1) + 1} \\
       &= \frac{726z - 725}{66z-65} \\
       &= \frac{726 - 725z^{-1}}{66-65z^{-1}}
\end{align}
$$

Which converted to a difference equation is:

$$
V_{in}[n] \approx 0.0909V_{out}[n] - 0.0895V_{out}[n-1] + 0.9986V_{in}[n-1]
$$

As far as whether retrieving the original signal is possible partly depends on whether the original signal did not in fact include significant energy above $f_s/2$.


The field is knows as Optical Character Recognition (OCR) and it's a broad field.

Since it's clearly an academic exercise, you can keep it simple. Start by image segmentation : find shapes where a shape is just a set of connected pixels. [1] Then compare the shapes to known shapes. Don't expect an exact match, do an approximate match. What overall properties can you reliably determine for the shape? Aspect ratio? Holes in the middle?

[1] At this point, you can already draw circles around each shape that's not too big nor too small.


Optimal time for sampling symbol in digital receivers

You are mixing up two different notions of sampling.

In a digital
communications system, the received (analog, continuous-time) signal
is passed through an A/D converter so that the needed further
processing (e.g. matched filtering) can be implemented on a
programmable DSP processor or as a MATLAB or C++ program or
on a special purpose ASIC. The sampling that you show in your question
is the one being carried out in the A/D converter, and as long as
the sampling is above the Nyqvist rate, the fact that the samples
are not at the peaks of the signal is irrelevant: the original
signal is represented perfectly adequately by the samples regardless
of where the peaks in the signal are with respect to the sampling
instants.

The other notion of sampling has to do with the matched filter
output which is a sequence of digital data. Which of these
data should the decision-making device use to decide whether
(in the simplest case) whether a 0 or a 1 was transmitted?
The answer is what you quote: it is the datum when the
signal output of the matched filter has a peak. That is,
we are picking one value from the matched filter output
sequence and making the decision based on that.
Note that

We are not going to average the matched filter output (or
take a weighted average of the output) and make a decision on that
quantity.
We are not going to take the maximum value of the actual output
sequence (which may, because of noise, be at a different location
than the time when the signal output of the matched filter is
supposed to peak) and make a decision based on that.

For more information on matched filters (including explanations
of the above bulleted points), see
this other answer of mine.


Processes: Orthogonal, Uncorrelated, Statistically Independent

You got some definitions wrong. It's correct that orthogonality means that $E[XY]=0$. Uncorrelated means that $X-\mu_X$ and $Y-\mu_Y$ are orthogonal, i.e., $E[(X-\mu_X)(Y-\mu_Y)]=0$. If you work that out you should arrive at the equivalent condition $E[XY]=\mu_X\mu_Y$ for uncorrelatedness (not for independence!). Consequently, if at least one of the two RVs $X$ and $Y$ have a zero mean, then orthogonality implies uncorrelatedness and vice versa.

Statistical independence means that the joint PDF of two random variables can be written as the product of the individual PDFs:

$$f_{XY}(x,y)=f_X(x)f_Y(y)\tag{1}$$

Independence implies uncorrelatedness, but the opposite is generally not true. If $X$ and $Y$ are jointly Gaussian, then independence and uncorrelatedness are equivalent. Consequently, in the special case that $X$ and $Y$ are jointly Gaussian and at least one of them has a mean of zero, then orthogonality, uncorrelatedness, and independence are all equivalent.


Real-time resampling pitch-shifting (without corrective time-stretching)?

okay, the answer to the question is "no", from my understanding of what "real-time" in DSP means.  a real-time DSP process is one that can run forever without either starvation of input samples or overflow of input samples.  for every second of sound input, there is a second of sound output.  just speeding up the sound or slowing it down, to change the pitch, is an operation you can do on a segment of sound or a sound file.  but it is not a continuous operation that runs real time.  this is not because the CPU or DSP cannot handle the number of operations, it's because, either the input or the output will have no time when the other still has time remaining.

now, there are DAW functions that play back sound files faster or slower that recorded without changing pitch.  what they are really doing is playing back the sound file faster or slower (with an increased or decreased number of input samples) and the DSP algorithm is doing real-time and complementary resampling and time-scaling to get to the output.


Difference Between Convolution and Multiplication

Multiplication in your first sentence is term-by-term multiplication: $z[n] = x[n]y[n]$ for all $n$. Convolution, for discrete-time sequences, is equivalent
to polynomial multiplication which is not the same as the term-by-term
multiplication. Convolution also requires a lot more calculation: typically $N^2$
multiplications for sequences of length $N$ instead of the $N$ multiplications
of the term-by-term multiplication.

The key point of Fourier analysis is that term-by-term multiplication in one
domain is the same as convolution in the other domain.  So, in order
to calculate the results of a convolution, you can either do it directly,
using $N^2$ multiplications, or transform to the other domain, do a term-by-term
multiplication, and transform back. This requires two transformations
to go from one domain to the other, $N$ multiplications in the other domain,
and one inverse transformation to come back to the domain
where the convolution result is needed. This more complicated process
can, in fact, require less computation because the transformations can
be done very efficiently via the Fast Fourier Transform (FFT) algorithm
which requires about $N \log N$ multiplications. So, compare the computational
effort in calculating the three transforms and doing the term-by-term
multiplication
to $N^2$ to see if the FFT gives you a more efficient method of computing
a convolution.


If your signal's in-phase and quadrature components are $x_I(t)$ and $x_Q(t)$, respectively, then your quadrature mixed signal is

$$s(t)=x_I(t)\cos(\omega_ct)-x_Q(t)\sin(\omega_ct)\tag{1}$$

where $\omega_c$ is the carrier frequency in radians per second. You can also use a '$+$' sign in $(1)$, that's just a matter of convention.


Combining two wave generation formulas does not work as expected (MATLAB/OCTAVE)

OK, got this solved by @omersayli (dsprelated).
I had put the variable t to a wrong location in both (triangle and square wave) formulas. Here's the correct formula for calculation of saw:
saw = -(2.*m.*asin((s-1)*cos((pi*f0*t)/w)).*acot(s.*csc((pi*f0*t)/w)))./pi^2.-sign(m).+m./2;


kernels to compute second order derivative of digital image

Convolution $(*)$ is associative, so

$I_{xx} = K_{sobel} * (K_{sobel} * I) = (K_{sobel} * K_{sobel}) * I$

meaning the convolution of a 1st order Sobel filter kernel with another 1st order Sobel filter gives the filter kernel for a second order filter. So e.g. a 2nd order horizontal Sobel filter kernel is:

$\left(
\begin{array}{ccc}
 -1 & 0 & 1 \\
 -2 & 0 & 2 \\
 -1 & 0 & 1 \\
\end{array}
\right) * \left(
\begin{array}{ccc}
 -1 & 0 & 1 \\
 -2 & 0 & 2 \\
 -1 & 0 & 1 \\
\end{array}
\right)=\left(
\begin{array}{ccccc}
 1 & 0 & -2 & 0 & 1 \\
 4 & 0 & -8 & 0 & 4 \\
 6 & 0 & -12 & 0 & 6 \\
 4 & 0 & -8 & 0 & 4 \\
 1 & 0 & -2 & 0 & 1 \\
\end{array}
\right)$

And for $I_{xy}$:

$\left(
\begin{array}{ccc}
 -1 & 0 & 1 \\
 -2 & 0 & 2 \\
 -1 & 0 & 1 \\
\end{array}
\right) * \left(
\begin{array}{ccc}
 -1 & -2 & -1 \\
 0 & 0 & 0 \\
 1 & 2 & 1 \\
\end{array}
\right)=\left(
\begin{array}{ccccc}
 1 & 2 & 0 & -2 & -1 \\
 2 & 4 & 0 & -4 & -2 \\
 0 & 0 & 0 & 0 & 0 \\
 -2 & -4 & 0 & 4 & 2 \\
 -1 & -2 & 0 & 2 & 1 \\
\end{array}
\right)$

If you want a 3x3 kernel, you can use the 1st order derivative filter $\left(
\begin{array}{cc}
 -1 & 1 \\
 -1 & 1 \\
\end{array}
\right)$, which leads to the second order filters:

$\left(
\begin{array}{cc}
 -1 & 1 \\
 -1 & 1 \\
\end{array}
\right)*\left(
\begin{array}{cc}
 -1 & 1 \\
 -1 & 1 \\
\end{array}
\right)=\left(
\begin{array}{ccc}
 1 & -2 & 1 \\
 2 & -4 & 2 \\
 1 & -2 & 1 \\
\end{array}
\right)$

and

$\left(
\begin{array}{cc}
 -1 & -1 \\
 1 & 1 \\
\end{array}
\right)*\left(
\begin{array}{cc}
 -1 & 1 \\
 -1 & 1 \\
\end{array}
\right)=\left(
\begin{array}{ccc}
 1 & 0 & -1 \\
 0 & 0 & 0 \\
 -1 & 0 & 1 \\
\end{array}
\right)$


I have been facing this problem of phase synchronization with my PCM/PM/SP-L signals as well, and found a temp solution by filtering away the residual carrier.

For carrier offset correction, there are two gnuradio blocks that work, the FLL Band Edge and PLL Carrier Tracking, but I found that the FLL Band Edge blocks works better together with the built-in Costas Loop block when dealing with PM signals with a residual carrier. This is completely from anecdotal experience though, and I have not done any proper characterization of the blocks performance.

One hacky solution of the phase synchronization problem might be to filter away the carrier of the signal before feeding it into the Costas Loop. This filter step should be after carrier offset correction, before Costas Loop. Since the presence of the residual carrier makes the IQ constellation of a Binary PM signal offset from the origin, filtering away the carrier will center the PM signal about the origin, and the Costas Loop PLL can lock its phase to that of a standard BPSK.

Rough Flow Graph:

PCM/PM RX -> FLL Band Edge -> Highpass Filter -> PFB Clock Sync -> Costas Loop -> BPSK Demodulator


Is there any new scientific capability to build non-causal filters in real world?

is there any new scientific theory which could yield to making none causal FIR Filters in the real world?

No. If anything, scientific theories seem to be challenging determinism itself, which is what could allow you to peer ahead in the future. The key here is the last phrase of that quote ("...many times the signal is stored and then processed...").

If you already had acces to the true $x[n+1]$ at the level of observation of a signal, why not use a causal filter anyway? What is so special about $x[n+1]$?

For instance, take a very simple finite difference filter. You can definitely express it as:

$$ y[n] = x[n+1] - x[n]$$

But you would not be able to apply it in real time. The objective here is to take the difference between two samples. If you wait for $x[n+1]$ to happen, you can simply:

$$ y[n] = x[n] - x[n-1]$$

If it was possible for this difference to work in the presence of $x[n+k]$, then $y$ could start producing some output before $x$ even begins to occur.

You can observe exactly what this would sound like in a time reversed reverb. Here is an example of that.

To try this yourself, load a sound in Audacity, apply "Reverse", "Reverb" and finally "Reverse" again, before hitting "Play".

In reverberation, present samples bear some resemblance to past samples. If you stand at the far end of a large room, sounds from sources with a direct line-of-sight to your ears reach you much faster than sounds that have to bounce off of the walls. The sounds that bounce off of the walls travel a longer distance and therefore reach you later. This is what is causing the time difference of different reflections in reverberation.

By reversing the sound clip, you make the current sample to bear resemblance to its future values, thus starting to hear reflections of an event that has not yet happen.

This is only possible with signals that are stored.

Hope this helps.


How to Pick a Journal for Publishing a Research Work?

Getting work published in a high impact factor journal is quite a slog. If you're starting from not having presented at a conference etc. before, that's going to be an even bigger ask.

However, it's not impossible.

For IEEE journals, short notes about signal processing topics can be sent to the IEEE Signal Processing Letters journal.

However, starting there might be too much.

What I'd suggest is to look at starting with an online journal such as arxiv or PLOSone.

Even before submitting there, I'd suggest posting a question here that might be relevant to the situations where the algorithm might be better than others, and also posting an answer with your algorithm.

The aim with any of these submissions is to get feedback about your work and improve it or its presentation.


How can I find SNR, PEAQ, and ODG values by comparing two audios?

I recommend calculating the correlation coefficient; the correlation coefficient will directly tell you how much B is (linearly) different from A, and you can compute the SNR in dB from the correlation coefficient $\rho$ using

$$SNR = 10log_{10}\left(\frac{\rho^2}{1-\rho^2}\right)$$

Given you have a noise-free version of your waveform, this would be a perfect application for the approach above. In this case "noise" is any deviation from your noise-free waveform once the two signals are properly scaled (which the computation will do).

See this similar post where I detail the correlation coefficient calculation:

Noise detection


Are you using scipy by any chance ?

If yes , this might help

https://scikit-dsp-comm.readthedocs.io/en/latest/_modules/sk_dsp_comm/digitalcom.html

check for functions "QPSK_tx"  and "QPSK_bb"


The phase of signal when it passes a frequency divider?

That can't generally be said. It depends on the actual design of the device.

In fact, it's not going to be constant for all frequencies, because that's physically pretty much impossible to achieve.

You'll find multiple descriptions in literature for three- and four-port devices like these, and you'll have to learn how to infer the relative phases for different frequencies from the measurements that the manufacturer of a divider will provide to you.


Segmentation : Berkeley Segmentation Dataset and Benchmark

Stereo : Middlebury Stereo Datasets

Shape from specular flow - Shape from specular flow dataset

Also, check out this link : http://clickdamage.com/sourcecode/cv_datasets.php


Math.Net and alglib returning different FFT outputs by default

I worked it out by myself, after reading a bunch of posts mentioning different methods of FFT output scaling. I still find this aspect of FFT processing heavily unsdocumented everywhere. I have not yet found any reliable source that explains what is the use of these scalings, which fields of sciences or what processing methods use them.

I have yet found out three different kinds of scalings, regarding the raw FFT output (complexes' magnitudes). This means multiplying them by:
1. 1/numSamples
2. 2/numSamples
3. 1/sqrt(numSamples)
4. (no scaling)

MathNet.IntegralTransforms.Fourier.Forward function (and according to various posts on the net, also possibly Matlab and Maple) by default, uses third one. Which results in the better distinguishable graphical output when using logarithmic colouring, in my opinion.

I would still be grateful if you know something and share your ideas, or if you can reference a good paper explaining on these.


Plotting Frequency Response Magnitude and Phase for first order all pass filter

As pointed out by Hilmar in a comment, you should remain in the discrete domain. This is most easily done by using the functions freqz() and zplane():

b = [.8,1]; a = [1,.8];
[H,w] = freqz(b,a,2048);
subplot(2,1,1), plot(w/pi,abs(H));
subplot(2,1,2), plot(w/pi,angle(H));
figure, zplane(b,a)

Note that for a first-order all-pass filter it's also quite straightforward to come up with analytical solutions for the magnitude and phase of the frequency response.


What happens if an harmonic in a band-limited signal at Nyquist frequency is added to a 90 degrees out-of-phase replica of it?

No.  Adding two waves at the same frequency results in another wave of the same frequency.  Depending on the phases, there can be anywhere from complete constructive interference to maximum destructive interference.  If the amplitudes are the same there can be complete cancellation.  This applies to the Nyquist frequency too.

See the math from my answer to FFT Frequency bin relationship.

Hope this helps.

Ced


Better understanding of downsampling (decimation) and upsampling (interpolation)

There are few problems in the code and question.

In your question 'given a sine wave at a given frequency for example, how would you proceed to downsample it to a lower frequency' this statement is not right. We are not reducing the frequency of sinewave, we are reducing the sampling frequency with which we sampled it.

Second, the code just inserts zeros alternatively. It doesn't downsample (means remove those samples altogether). This means, it is still sampled with original sampling frequency (8kHz) but now with alternate zeros. And what is the effect of alternate zeros in signal? It creates a replica of spectrum.
Also, in the code, the signal frequency is exactly half of sampling frequency, which causes sampling to happen at zero crossings. You can plot the signal $x$ and see.

In the code, I changed $f$ to 3000, and can see 2 peaks at 3000Hz and 5000 ($fs -  f$)Hz. For the downsampled signal, apart from original signal, there is 2 peaks corresponding to replicated spectrum.


Determining linearity, causality, memory, and time invariability from a picture of the signal

Basically you can't. In basic information theory a system has an input and output and the terms you are asking about are defined in terms of the relationship between output and input. It's difficult to determine what the input to saccadic eye movements would be. It's mostly an unconscious process controlled by the neural periphery.

In order to answer your question you would first define an appropriate input signal and than look at a screen shot that shows both input and output simultaneously.

FWIW you could potentially argue that it's causal simply because it is a real physical system and all physcial systems are causal.


Can use of Fourier transform be minimized completely with the help of Laplace and Z transform?

The short answer is yes, if you have the Laplace or Z-transform of a function you do not need the Fourier transform.

This is because the CFT is a special case of the Laplace transform and the DTFT is a special case of the Z transform. The Fourier transform is used to find the complex sinusoids that compose a function, whereas the Laplace transform finds all complex exponentials that compose a function.

Knowing this, it is simple to derive the Fourier transform from the Laplace transform. The Laplace transform is defined as:

$$
\int_{-\infty}^{\infty}e^{-s}f(t)\,dt\\
s = \alpha+j\omega
$$

To find the Fourier transform from the Laplace transform, all you must do is assert $\alpha = 0$ for all $\omega$ and $t$. This evaluates the Laplace transform along the imaginary axis only, as any component lying directly on the imaginary axis is a complex sinusoid. This leaves us with:
$$
\int_{-\infty}^{\infty}e^{-j\omega}f(t)\,dt\\
$$

Which is in fact the Fourier transform.

A similar relationship applies between the DTFT and the Z-transform. The Z-transform is defined as:

$$
\sum_{n=-\infty}^\infty x[n]Z^{-n}
$$
Where $Z$ is a complex variable. To find all complex sinusoids in the discrete world we must constrain $Z$ to fall on the unit circle:
$$
Z = e^{j\omega}
$$

Which leaves us with

$$
\sum_{n=-\infty}^\infty x[n]e^{-j\omega n}
$$
Or the discrete time Fourier transform.

The reason the Fourier transform exists is because there are a lot of applications that only need the real frequencies of a function. If only the real frequencies are needed there is no point calculating the rest of the Laplace domain. Why do more work than necessary?

For more information check out the Wikipedia articles https://en.wikipedia.org/wiki/Discrete-time_Fourier_transform#Relationship_to_the_Z-transform
and https://en.wikipedia.org/wiki/Laplace_transform#Fourier_transform And I would also strongly recommend the text book Linear Systems and Signals by B.P. Lathi for a much more indepth discussion on Laplace, Z, and Fourier transforms.


follow up to previous question " is this an impulse response or step response"

Let me provide an answer,

1-) In the answer I have used the notation $x(t)=\delta(t)$ and $x(t)=u(t)$ to indicate that the input to the system is a dirac impulse and unit-step function respectively. The correct notation should be $x(t) \equiv  \delta(t)$ which means identical, by definition. Triple equality sign is almost always omitted and the double equality sign is used instead in signal processing context. (it's therefore not asking for the solution of an equation, which is generally understood from double equality sign).

2-) The subject of differential and difference equations is rich. According to the usual literature, continuous time-systems are described by differential equations and discrete-time systems are described by difference equations. When a continuous-time signal is sampled and the associated continuous system is transformed into a discrete-time system the corresponding differental equation is also discretized and converted into an equivalent difference equation. One application of such discretization is the usual Kalman filering equations.

I've analysed the continuous-time leaky integrator model whose input-output relationship was assumed by definition as $y(t) = \mathcal{T} \{x(t) \} = \int_{-\infty}^{t} e^{-\alpha(t-\tau)} x(\tau) d\tau $. Given this definition of the system the corresponding continuous-time LCCODE is obtained as in the answer. In the web-link only the differential equation relation is provided as $ dy/dt = s - y(t)/k $ where $y(t)$ is the water level (the output) $s=s(t)$ is the input and $1/k$ is the leakage constant, the differential equation is then $$ \frac{dy}{dt} + \frac{1}{k} y = s $$ or equivalently $$ \frac{dy(t)}{dt} + \frac{1}{k} y(t) = s(t) $$ This differential equation is the same as $ y' + \alpha y = x$ after replacing $\alpha = 1/k$ and $s(t) \equiv x(t)$.

On the other hand in your original question, you refer to a discrete differential equation that I cannot clearly understant whether it's a difference equation or a differential equation (the terms discrete and differential used together as a mix) or it's just another application of discretization as I defined above.

However the following excerpt from the original web-reference indicates to me that:

  We can easily simulate this leaky integrator in discrete steps of time
  by changing the value of 'y' on each step according to the equation
  above:

There is only a continuous-time differential equation and it's solution is numerically computed by the software, which is indicated as the simulation of the ODE using discrete steppings (the numerical computation)

So as I understood given the continuous time differential equation
$$ y'(t) + \alpha y(t) = x(t) $$ one can reach the numerical version as
$$ \frac{d y}{dt} = -\alpha y(t) + x(t) $$
$$ dy = ( -\alpha y_t + x_t) dt $$

Where $y(t)$ is replaced by $y_t$ just as a matter of notation. So they are the same thing. The web-link numerically solves the contiuous differential equation for various input stimulus.


Spectral Repair is more complex than just muting a handful of harmonics.

But first, let's just clarify that as far as Linear Time Invariant (LTI) systems are concerned, their processing can be applied, equivalently, either in the time domain (via the operation of convolution) or the frequency domain via either the overlap-add or overlap-save methods.

iZotope's Spectral Repair performs interpolation. That is, to remove an unwanted portion of a recording, it will try to "mimic" the characteristics of the spectrum in either ends of the unwanted region.

Consider the following image:

Can't we just remove the railings by copying and pasting the grass from another (and very similar) part of the image and adjust its brightness? Couldn't we even replicate how the green colour varies in the grass portions and just generate some "new" grass for the regions of the railings?

Is it going to look odd? Possibly, if it was overdone, but for small regions it may be good enough to fool the eye.

This is what the "Replace" mode (and "Attenuate") does in iZotope, but instead of grass and railing there is background noise and local disturbance.

It's not so much attenuation (or "setting harmonics to zero") as "masking" or "hiding" the unwanted sound, since it attempts to make up a good enough patch of harmonics to bury it in by "looking at" the surroundings of the disturbance.

For more information on interpolation please see this link. For an example of "learning the profile and applying a filter" please see this link and this link.

Hope this helps.


What does the physical meaning of $\mathbf{U} \Sigma \ \mathbf{V}^{*} $ in the MIMO system?

Singular value decomposition (SVD) is a method of decomposing a matrix into two unitary matrices and a non-negative definite matrix of singular values.

In the context of MIMO communications, the $\mathbf{H}$ matrix defines the channel; it doesn't specify anything directly about the transmitter or the receiver.

This paper has the following diagram which shows one interpretation of the $\mathbf{U}$ and $\mathbf{V}$ matrices.  The SVD allows the channel to be decomposed into several parallel paths, with no crosstalk between them.


Is the common FIR Filter delay line (see pic) also valid on purely logical FPGA implementations?

I'm not sure if I understood our question right, and therefore sorry if I point something wrong, but the delay in the delay line comes from signal sampling and not from computational processing type. Therefore whether you realize that signal flow graph in a von Neumann type serial architecture or in aparallel FPGA based architecture doesn't matter. You will always have the delay per input sample as dictated by your systems sampling rate.

On the other hand of course a parallel implementation would take much less clock cycles, based on the amount of parallelisation, to compute per sample output campared to a serial architecture of a general purpose CPU, however as I said the delay in the delay line is nothing related with computational delay of the hardware system.

Finally the signal flow graph yields a pictorial representation of processing structure, you don't really implment those delays. They are represented by registers or memory locations, as possibly you know, which are updated either per sampling clock or via buffer copies which provides that delay effect represented in the flow graph.


Audio steganography using phase encoding technique

Here is the algorithm:

Divide an input signal $S$ into $n$ consecutive blocks $S_i$ ($i\in\{1,\ldots,n\}$) of length $l$ each
Compute magnitude $A_i$ and phase $\phi_i$ of each block using FFT where A(i)=abs(fft(S{i})) and Phi(i)=angle(fft(S{i})) in MATLAB.

For example, here are the first 3 elements of $\phi$:

(Recall that FFT of real signal has odd symmetry which is clearly visible on plots)
Compute phase differences $\Delta\phi_i=\phi_i-\phi_{i-1}$ for $i\in\{2,\ldots,n\}$
To encode your binary data $d$ of length $m:m<l$ assign
$$\phi_{data}[i]=
\begin{cases}
\pi/2,& d_i=0\\
-\pi/2,& d_i=1
\end{cases}$$ for $i\in\{1,\ldots,m\}$.
Replace following elements of first phase sequence, $ {\phi_1} $,  with ${\phi _{data}}[i]$ for $i\in\{1,\ldots,m\}$ as:
$${\phi '_1}[ {L/2 - m + i}] = {\phi _{data}}[i]\,$$
To maintain odd symmetry property of DFT, repeat same progress as:
$${\phi '_1}[{L/2 + 1 + i}] =  - {\phi _{data}}[{m + 1 - i}]$$
To maintain phase differences sequentially reassign $\phi '_i=\phi_{i-1}+\Delta\phi_i$ for $i\in\{2,\ldots,n\}$
Here are the first 3 elements of $\phi '$ after this step:

(Note how encoded $d$ has affected next blocks)
Reconstruct the signal using inverse FFT applied to each block $A_i exp(j \phi '_i)$ (where $j$ is imaginary unit) and joining all blocks together.

Explanations to many details like the length of blocks and why particular steps are taken are given in the documents you and endolith linked to.

A minimal IPython Notebook version: github.com/danylo-dubinin/secret_in_wav.


How to sketch output PDF given the transformation?

There is no impulse because there is zero probability mass at $0$. As the problem itself shows, it is perfectly all right for a pdf to make sudden jumps in value; note that the graph of $f_X$ (yes, I do mean $X$, not  $Y$) jumps in value at $0$ and the value of $f_X(0)$ cannot be determined from the graph. $f_X(0)$ could have value $\frac 14$ or $\frac 12$ or any value in between (or any nonnegative value for that matter)


What is the frequency response function for an input through a differentiator and mixer?

The output of the differentiator is $Y(\omega)= j\omega X(\omega).$

Multiplying by the sine in the time domain is convolving with the Fourier Transform of sine in the frequency domain:

$$\mathcal F\left\{\sin(w_ct)\right\} = \frac{1}{2j}(\delta(\omega+\omega_c)-\delta(\omega-\omega_c))$$

Which results in the following expression for the output of the mixer, as a function of $Y(\omega)$:

$$\frac{1}{2j}(Y(\omega+\omega _c)- Y(\omega-\omega _c))$$

substituting back for $Y(\omega)$ in terms of $X(\omega)$:

$$\frac{1}{2j}(j(\omega+\omega _c) X(\omega+\omega _c)- j(\omega-\omega _c) X(\omega-\omega _c))$$

which is

$$\frac{1}{2}(\omega+\omega _c) X(\omega+\omega _c)- (\omega-\omega _c) X(\omega-\omega _c))$$


How to validate that an audio algorithm result is independent of microphone device?

There isn't a universal pre-processing chain, but here are some common ones that I'd look out for:
1) Noise Gating
If a signal is below a certain threshold, the mic is muted. To identify: You'll see a VERY quiet noise floor, followed by a sudden spike in level.
2) Automatic Gain Control (AGC)
Mic volume modulates based on the signal level. This will be difficult to identify unless you have multiple recordings from the same system, which will show modulation in the noise floor based on the signal level.
3) Compression
This one should be self-evident. You'll see a very steady signal level that makes the waveform look like a "block", i.e. not much changes in level throughout the entire recording.
These are all relatively simple real-time processing examples. There may also be more sophisticated post-processing applied in some cases, such as noise suppression.


Impulse Invariant method for digital filter design

This is just as it turns out when you do the math. The discrete-time Fourier transform (DTFT) of the sampled continuous-time impulse response $h(t)$ is
$$H_d(e^{j\omega T})=\sum_nh(nT)e^{-jn\omega T}\tag{1}$$
With
$$h(nT)e^{-jn\omega T}=\int_{-\infty}^{\infty}h(t)e^{-j\omega t}\delta(t-nT)dt\tag{2}$$
this can be written as
$$\begin{align}H_d(e^{j\omega T})&=\sum_n\int_{-\infty}^{\infty}h(t)e^{-j\omega t}\delta(t-nT)dt\\&=\int_{-\infty}^{\infty}\left[h(t)\sum_n\delta(t-nT)\right]e^{-j\omega t}dt\\&=\mathcal{F}\left\{h(t)\sum_n\delta(t-nT)\right\}\\&=\frac{1}{2\pi}H(\omega)\star\frac{2\pi}{T}\sum_k\delta\left(\omega-\frac{2\pi k}{T}\right)\\&=\frac{1}{T}\sum_kH\left(\omega-\frac{2\pi k}{T}\right)\tag{3}\end{align}$$
where $H(\omega)$ is the Fourier transform of $h(t)$, and $\star$ denotes convolution. From $(3)$ we see that the DTFT of the sampled impulse response equals the sum of shifted spectra of $h(t)$, scaled by $1/T$.
If we assume that $H(\omega)$ is approximately band-limited and that $T$ is chosen sufficiently small such that aliasing becomes negligible, we obtain the approximation
$$H_d(e^{j\omega T})\approx\frac{1}{T}H(\omega),\qquad |\omega|<\frac{\pi}{T}\tag{4}$$
For the step-invariance method, we use samples of the step response instead of samples of the impulse response, and we obtain a relation analogous to $(3)$ between the DTFT $G_d(e^{j\omega T})$ of the step response of the discrete-time system, and the Fourier transform $G(\omega)$ of the continuous-time step response:
$$G_d(e^{j\omega T})=\frac{1}{T}\sum_kG\left(\omega-\frac{2\pi k}{T}\right)\tag{5}$$
In order to obtain the frequency response $H_d(e^{j\omega T})$ we multiply $(5)$ by $1-e^{-j\omega T}$, because the impulse response is obtained by computing a first-order difference of the step response:
$$H_d(e^{j\omega T})=\left(1-e^{-j\omega T}\right)G_d(e^{j\omega T})=\frac{1-e^{-j\omega T}}{T}\sum_kG\left(\omega-\frac{2\pi k}{T}\right)\tag{6}$$
For frequencies that are small compared to the sampling frequency, i.e., for $|\omega T|\ll 1$ we obtain from $(6)$
$$\begin{align}H_d(e^{j\omega T})&\approx\frac{1-(1-j\omega T)}{T}\sum_kG\left(\omega-\frac{2\pi k}{T}\right)\\&=j\omega \sum_kG\left(\omega-\frac{2\pi k}{T}\right)\tag{7}\end{align}$$
If we again assume that aliasing can be neglected, we arrive at
$$H_d(e^{j\omega T})\approx j\omega G(\omega)=H(\omega),\qquad |\omega|<\frac{\pi}{T}\tag{8}$$
From $(8)$ we see that, unlike for the impulse invariance method, the step invariance method doesn't involve a scaling of the continuous-time frequency response.


Split and Merge image segmentation algorithm in OpenCV

I asked Google again for you, but I did manage to find some hits in the end. There is already a very good question on stackoverflow concerning the exact same thing you are interested in.

There is a very nice explanation of split-and-merge provided in one of the answers, as well as simplified pseudocode.

The other answer provides a link to the implementation (although I did not check it, I'm not sure which language the implementation is in).

If you want to use split-and-merge as a reference segmentation technique to compare your idea with, other segmentation methods offered by OpenCV might also be interesting:

functions useful in image segmentation from the OpenCV manual
a chapter dealing with image segmentation in OpenCV from Gary Bradski; Adrian Kaehler: Learning OpenCV
article about Pyramid segmentation
Mean-shift segmentation: stackoverflow answer, scientific article
Watershed segmentation seems fairly well explained by Wikipedia

In the end, you can choose one of the already offered methods as the reference method after all, or, if you want to present amazing, all-encompassing comparisson, you might also want to use more than one reference method


Causality and ROC of a stable LTI system

This looks like homework, so I'll give you a hint to help you figure this out by yourself.

There are 2 poles: one at $z=3$ (i.e., outside the unit circle) and one at $z=\frac13$ (i.e., inside the unit circle). For such a system there are three possible ROCs; one corresponds to an anti-causal system, another one to a causal system, and a third one to a non-causal system with a two-sided impulse response. But only one of those ROCs includes the unit circle, and that's the ROC corresponding to a stable system.


What is the impulse response of SIMO channel and its equivalent in oversampling?

Fractional sampling and SIMO are not exactly similar, SIMO system uses multi anntenna at receiver but Fractional sampling can be used with one antenna at receiver.  It's said that, Fractional sampling is often used to convert the SISO system into SIMO system without using Multi Antennas at receiver. How?

Suppose a SISO system is used, where you use single Tx antenna with single Rx antenna, Fractional sampling means you oversample the received signal by G, where G is  Oversample factor, that means your sampling rate becomes Ts/G, and Ts is the sampling rate used at transmitter. So, that means you don't change the structure of transmitter, but you oversample the receiver. How?

Suppose the received signal is :

where y,u,h are received signal, transmitted signal and channel respectively.

If y(t) is sampled with a rate of Ts/G instead of Ts, where G is the oversampling ratio, and "Ts is the sampling rate which was used at transmitter:  its polyphase components can be expressed as:

where g = 1, ..., G.

each polyphase component should be demodulated in a separated demodulator (This the one of disadvantage of using Fractional sampling). So instead of having one demodulator, we will have multi demodulators.

Coming back to the previous equation,

and the channel hg[n] is a matrix as below:

which means,

So now:

It seems similar,  but I can't confirm or give more details, I am not sure, AlexTP maybe can explain it better than me.

thank you and sorry for the notations, because I don't know how I can write the mathematical equations here.


Check a music production software (commercial: Logic/ Ableton/ Pro Tools. Open Source: Ardour, Mu Lab) or otherwise audacity. Add a stereo widening tool plugin.
For that check KVR Audio
stereo/ imaging/ widening
stereo imaging more close search


A system is causal if its output depends only on the current input and past inputs (and not on future inputs).  As a consequence, if the system is a linear, time-invariant (LTI) system (which input/output relationship can be completely characterized by an impulse response), that impulse response, $h(t)$ must be zero for all time $t \lt 0$.

Some people define a causal signal, $x(t)$, to be one that can be the impulse response of a causal system: it is zero for all time $t \lt 0$.


Decomposing an Image into High Frequency and Low Frequency Parts Using Bilateral Filter

Since the Bilateral Filter is Spatial Invariant and Non Linear (The output is linear combination of the data, yet the weights depends on the data in Non Linear fashion) it can't be formulated as a Spatially Invariant Filter / Convolution.

So the Low Pass / High Pass concept in it classic form doesn't hold here.

Yet you can still employ the Image Pyramid concept here with your operator being the Bilateral filter and everything will work perfectly (Namely perfect reconstruction).

It is actually very effective.


What about just integrating the energy in the high frequencies?
You can do that after applying High Pass filter which is equivalent of weighting high frequencies differently.


Basics of leakage phenomena in DFT and its approximation with sinc function

No, the two peaks don't need to be the same. The discrete-time Fourier transform (DTFT) (samples of which are computed by the DFT) is also not symmetric with respect to the sinusoid's frequency (as explained in Robert's comment). If you compute the DTFT of a truncated sinusoid

$$x[n]=\sin(\omega_0n)(u[n]-u[n-N])$$

you get

$$X(e^{j\omega})=\frac{1}{2j}\left(e^{-j(\omega-\omega_0)\frac{N-1}{2}}\frac{\sin\left(\frac{(\omega-\omega_0)N}{2}\right)}{\sin\left(\frac{\omega-\omega_0}{2}\right)}-e^{-j(\omega+\omega_0)\frac{N-1}{2}}\frac{\sin\left(\frac{(\omega+\omega_0)N}{2}\right)}{\sin\left(\frac{\omega+\omega_0}{2}\right)}\right)$$

The first term in the parentheses has its main contribution around $\omega=\omega_0$, whereas the second term has it around $\omega=-\omega_0$. Note that the two functions are no sinc functions; only for $\omega$ very close to $\omega_0$ (or $-\omega_0$) can they be approximated by a sinc function. And that still wouldn't take into account their overlap.

So your results are correct, and the plot below shows your DFT on top of the DTFT as computed above:

Here's the Matlab/Octave code for the DTFT plot:

N = 32;
k = 8.5;
fs = 32000;
w = 2 * pi * linspace(0,1,500);
w0 = 2*pi*k/N;
X = exp(-1i*(w-w0)*(N-1)/2) .* sin(N/2*(w-w0))./sin((w-w0)/2);
X = X - exp(-1i*(w+w0)*(N-1)/2) .* sin(N/2*(w+w0))./sin((w+w0)/2);
X = X/(2*1i);
plot(w*fs/(2*pi),abs(X))


Consider looking at this StackOverflow answer which provides the full code for creating a Butterworth bandpass filter. In the sample code the answerer provides, the filter is applied to a manually constructed simplistic signal x. If you substitute this signal with your own NumPy array, it should achieve the intended effect.

WAV files can be quickly made into NumPy arrays, consider using the librosa module. librosa.core.load will be what you use. To write out the filtered signal, simply use librosa.output.write_wav.


Question about convolution of shifted function-simple

Considering an LTI continuous-time system with impulse response $h(t)$, its response to any valid input $x(t)$ can be computed with a convolution integral as:
$$y(t) = x(t) \star h(t) = \int_{-\infty}^{\infty} x(\tau) h(t-\tau) d\tau$$ and equivalently
$$y(t) = h(t) \star x(t) = \int_{-\infty}^{\infty} h(\tau) x(t-\tau) d\tau$$
due to the commutativity of the convolution operator.
Hence, shifting the output $y(t)$ by $k$ could be shown as
$$y(t-k) = \int_{-\infty}^{\infty} x(\tau) h(t-k-\tau) d\tau$$ or equivalently as
$$y(t-k) = \int_{-\infty}^{\infty} h(\tau) x(t-k-\tau) d\tau$$
Where, the first one equals $y(t-k) = x(t) \star h(t-k)$ and the second one equals $y(t-k) = h(t) \star x(t-k)$.


Changes in frequency are equivalent to a linear phase sweep, so you're seeing exactly what you'd expect. Consider a reference sin wave with frequency $\omega_1$, and a signal with frequency $\omega_2=\omega_1+\Delta\omega$. Well,
\begin{align}
\sin\omega_2 t &= \sin((\omega_1 + \Delta\omega)t) \\
 &= \sin(\omega_1t+\Delta\omega t) \\
 &= \sin(\omega_1 t+\phi(t))
\end{align}
Where $\phi(t) = \Delta\omega t$. This is the line you're seeing in the unwrapped phase.

So if you do lock in at $\omega_1$, your're going to pick up that frequency component, but the frequency offset will be seen in the phase varying linearly with time.


Is system $h(t) = u(t)-u(t-1)$ LTI or not?

If a system is described by an impulse response $h(t)$, then by definition it is an LTI system. Because impulse response of the form $h(t) = \mathcal{T} \{ \delta(t) \} $ only exists for LTI systems.


What is the "pitch marker" in the LP-PSOLA algorithm?

The pitch marker indicates the beginning of each cycle of the waveform - or from a more physiological point of view the point at each the periodical motion of the glottis causes a sudden change in air flow.

A stable sound (in pitch and phoneme) will still have a marker every cycle - so these markers do not indicate any change.


Proof of transfer function factorization $\frac{b_0}{a_0} \frac{\prod_{k=1}^M (1-c_kz^{-1})}{\prod_{k=1}^N(1-d_kz^{-1})}$

It's a mere application of the fundamental theorem of algebra of polynomials in 1D, which states that any polynomial $p(x)$ of $x$ of degree $n$, can be written (expressed) in either increasing (or decreasing) powers of $x$ as in $$p(x)=a_0x^n+a_1x^{n-1}+...+a_{n-1}x+a_n = \sum_{k=0}^{n}{a_kx^{n-k}}$$  or as a product of its $n$ roots in the factors of $(x-r_k)$ as in $$p(x) = K(x-r_1)(x-r_2)...(x-r_n)=K\prod_{k=1}^{n}(x-r_n)$$ where each such factor represents a root of the polynomial $p(x)$.

The $n^{th}$ degree poynomial is assumed to have $n$ roots, either distinct or repeating or in complex and/or complex conjugate forms, depending on the coefficients $a_k$.

Applying this to the $Z$-Transform of the form discussed:
$$X(z)=\sum_{n=0}^{M}{a_nz^{-n}}=a_0+a_1z^{-1}+...+a_Mz^{-M}$$ this can be seen either as a polynomial in the variable $z^{-1}$, or after factoring out $a_0z^{-M}$ ($a_0$ being nonzero) it will be a polynomials in $z$ as given $$X(z)=a_0z^{-M}\left(z^{M}+\frac{a_1}{a_0}z^{M-1}+...+\frac{a_M}{a_0}\right)$$ And the paranthesis can be written in product form as:
$$X(z)= a_0z^{-M}\prod_{k=1}^{M}{(z-r_k)}$$
You can distribute the factor $z^{-M}$ into each product form as in: $$z^{-1}(z-p_m) = (1-p_mz^{-1}) $$ which yields the final required form as:
$$X(z) = a_0 \prod_{k=1}^{M}(1 - r_k z^{-1})$$ where each $r_k$ is a root.

When $X(z)$ is expressed as a rational traform as in $$X(z)= \frac{P(z)}{Q(z)} = \frac{\sum_{k=0}^{M}{b_kz^{-k}}}{\sum_{k=0}^{N}{a_kz^{-k}}} $$ apply this factorization to both of numerator and denominator to get the final form:$$X(z)= \frac{\sum_{k=0}^{M}{b_kz^{-k}}}{\sum_{k=0}^{N}{a_kz^{-k}}}=\frac{b_0}{a_0} \frac{\prod_{k=1}^{M}{(1-c_kz^{-1})}}{\prod_{k=1}^{N}{(1-d_kz^{-1})}}$$


Straight lines with an edge detection algorithm

You're probably looking for the Hough transform or one of it's extensions.

The simplest version of this transform is linear and appropriate for detecting straight lines.

In the transformed space (Hough space), angles and distances are found as points where curves intersect.

Libraries for calculating the Hough transform exist in

C++ - OpenCV (Has Python bindings)
Matlab - hough in the image processing toolbox
Python - Scikits-image

Related Questions and Answers

See this previous answer for more help on understanding the Hough transform.


Read the original paper: Schmidt, R. O. "Multiple Emitter Location and Signal Parameter Estimation." IEEE Transactions on Antennas and Propagation. Vol. AP-34, March, 1986, pp. 276–280

You may also want to look up "Pisarenko's Method", "Prony's Method" and read about related problems such as ESPRIT (Roy, R.; Kailath, T. (1989). "Esprit - Estimation Of Signal Parameters Via Rotational Invariance Techniques". IEEE Trans. Acoustics, Speech and Signal Processing 37 (7): 984–995.)

You can find a lot of this in statistical signal processing books like "Statistical Digital Signal Processing and Modeling" by Monson Hayes and look at the cited references therein.  Another reference is Mathematical Methods for Signal Processing by Moon & Stirling.


Mismatched offset when computing cross-correlation using FFT

You're missing an fftshift to get the correct alignment:

plot(lags,fftshift(result));


Spectrogram from scipy.signal with Python

The input to your spectrogram is too short in terms of time.  To get 2 Hz plot resolution, you need an input vector long enough to represent at least 0.5 seconds (preferably 3X more) worth of samples at your sample rate.


Is the decimation process time invariant?

Sample rate conversion systems (expansion or decimation) are time-varying operators.

For your example system of decimation by $M$ :

$$ y[n] = T\{x[n]\} = x[Mn] $$

you can easily see that results of shift in the input and output are not the same; i.e.,

$$ y_1[n] = T\{ x[n-d] \} = x[Mn - d] $$

and

$$ y_d[n] = y[n-d] = x[M(n - d)] \neq y_1[n]$$


What does $[-\pi \; \pi]$ band-limited signal mean?

Now first of all concept of bandlimited signal exists in the context of analog signals not digital ones. Digital signals are always (by their definition and existance) bandlimited.

Second, frequency in analog signals is expressed either as $f$ in Hertz (cycle per second, repetation per second) or as $\Omega$ in radians per second, where they are related as $$\Omega = 2 \pi f$$ So for example a frequency of $f=1$ kHz would equivalently make a radian-frequency of $\Omega=2 \pi 1k \approx 6283$ rad/sec.

When a bandlimited analog signal $x_a(t)$ is sampled by a period of $T_s$ (and digitized afterwards in practice), the resulting discrete-time sequence $x[n]$ will have a discrete Fourier transform expression of
$$ H(e^{j\omega}) = \frac{1}{T_s} \sum_k X_a( \frac{\omega + 2\pi k}{T_s})$$

such that it will be frequency normalized to the discrete-time frequency range of $-\pi \leq \omega < \pi$. Where $\omega$ is the discrete-time frequency in radians per sample. Where $X_a$ is the continuous-time Fourier transform (frequency content) of the analog signal $x_a(t)$.

for $k=0$ the base period of $H(e^{j \omega})$ becomes $$ H(e^{j \omega}) = \frac{1}{T_s} X_a( \frac{\omega}{T_s}) $$ where the frequency normalization becomes more apparent as setting $ \Omega = \frac{\omega}{T_s} $ yields the those analog radian-frequencies for $X_a(\Omega)$ when $\omega$ of $H(e^{j \omega})$ is selected from $\omega \in [-\pi,\pi)$


Running an FFT filter on a large data set

You need to know the length of the impulse response of your bandpass filter, and use an FFT longer than 256 by that length (to prevent circular convolution boundary artifacts).  Combine the longer FFT results back into a sequence of length 256 windows by overlap-add or overlap-save.


Why does this manual bilinear transformation yield different results from Matlab's?

A couple things:

Before making the substitution $s = \frac{2}{T} \frac{z-1}{z+1}$, you need to prewarp the cutoff frequency by making the substitution:

$$
\omega_{c,w} = \frac{2}{T} \tan (\omega_c\frac{T}{2})
$$

where $\omega_{c,w}$ is the warped cutoff frequency. This is necessary because the bilinear transform maps the left-half plane in the Laplace domain (used in analog filter design) to the unit circle in the $z$-domain in a nonlinear fashion. Therefore, as you approach the Nyquist rate (digital frequencies of $\pm \pi$), the approximation to the analog filter prototype becomes inaccurate.

Also, the second parameter that you are passing to the butter function is the normalized cutoff frequency, not the sampling interval $T$. The normalized frequency used by that function is in the interval $(0,1)$ and is equal to the ratio of the desired cutoff frequency to the Nyquist rate:

$$
\omega_n = \frac{\omega_c}{2 \pi \frac{f_s}{2}}
$$

$$
\omega_n = \frac{\omega_c}{\pi f_s}
$$

$$
\omega_n = \frac{\omega_c T}{\pi}
$$


From your statement, i understand a discrete-time downsampling of a sequence $x[n]$ already converted to digital, without aliasing, previously.

So to prevent any aliasing on the downsampled sequence $y[n]$, you should first apply an anti-aliasing lowpass filter to the input sequence $x[n]$. The cutoff frequency of the lowpass filter is given by the downsampling ratio $D$ and is:

$$ \omega_c = \frac{ \pi }{D }$$ radians per sample. The gain of the lowpass filter is one.

The following Matlab\Octave code designs a linear phase, type-I, FIR lowpass filter of odd length $L = 2K+1$ with  a group delay of $K$ samples :

D = 8;     % downsampling factor
K = 256;   % group delay of FIR linear phase filter.
h = fir1(2*K, 1/D); % FIR, LPF, Type-I, L = 2K+1, wc = pi/D


Difference Between Linear Convolution and Circular Convolution

The difference applies only to the borders of the image.

In the linear convolution you assume the values of pixels beyond the border (examples being mirror of the image pixels, or 50% grey).

In the circular convolution (or DFT, product, IDFT), the pixels beyond the border are the pixels on the other side of the image, just as if you had a repeated tiling of the image.


why always structured Images are considered as 2D signals?

A 1D signal that is only noise with no structured feature to be seen is of no interest.

A 2D signal that is only noise or a superposition of uncorrelated 1D signal rows, with no 2D structured feature to be seen, is of no interest.


You can do it like this. I think the code speaks for itself:

x = [1, 2, 3, 4, 5, 6, 7, 8];   % some real-valued signal
n = length(x);                  % must be even!
n2 = n/2;                       % assume n is even
z = x(1:2:n)+j*x(2:2:n);        % complex signal of length n/2
Z = fft(z);
Ze = .5*( Z + conj([Z(1),Z(n2:-1:2)]) );        % even part
Zo = -.5*j*( Z - conj([Z(1),Z(n2:-1:2)]) );     % odd part
X = [Ze,Ze(1)] + exp(-j*2*pi/n*(0:n2)).*[Zo,Zo(1)];  % combine
X2 = fft(x);
max(abs(X-X2(1:n2+1)))  % 2.4492e-15


Why is no phase information lost when a received signal goes through a mixer?

The receiver's mixer introduces a phase-shift to the in-phase and quadrature components. In very general terms, assume we have a signal $$a(t)e^{j(2\pi f_0 t + \phi(t))},$$ which is mixed with a local oscillator of frequency $f_{lo}$ and phase $\delta$: $$e^{j(2\pi f_{lo}t + \delta)}.$$ The output signal is

\begin{align}
a(t)e^{j(2\pi f_0 t + \phi(t))}e^{j(2\pi f_{lo}t + \delta)}
&= a(t)e^{j(2\pi (f_0+f_{lo})t + \phi(t) + \delta)} \\
&= a(t)e^{j(2\pi (f_0+f_{lo})t + \phi(t))} e^{j\delta}.
\end{align}

In other words, even though the frequency has been shifted and a phase-shift equal to $\delta$ has been introduced, $a(t)$ and $\phi(t)$ are still present in the signal and that means that all the information carried by the original signal is still there.

Usually the phase-shift is corrected in the digital back-end in the receiver, or you can use a differential encoding (such as DBSK), in which the phase-shift becomes irrelevant.


Gradient echo MRI. Do the gradients have to be the same duration?

The question is what optimal is. It is correct that very often certain gradients do have the same duration in sequence diagrams - and actual imaging sequences (e.g. slice rewinder, read-out prewinder, and phase-encoding gradient). Normally, you would compute the longest duration for these gradients and recalculate the others to be of the same length by stretching them in time and reducing their amplitude. This may have benefits, e.g. the peripheral nerve stimulation may be reduced, and also the noise that the MRI creates. Also it reduces the stress for the gradient coils, i.e. heating and vibration. Certainly it is not always the best idea to stretch all gradients: When the fall below a certain amplitude, there may be inaccuracies due to system imperfections. Also the parasitic diffusion weighting of the gradients may be an issue in certain applications.

So the main reasons for certain design choices are:

reduction of patient discomfort
reduction of stress for the MR-system
reduction of image artifacts


Simulating noise in computed tomography reconstruction

The easiest, and probably most straight forward way, seems to add the noise in the measurement domain, hence the sinogram.


How to equalize an ISI channel when the transmitted symbols are unknown?

The classic approach is to include a "training sequence" along with the data. The training sequence is known to both transmitter and receiver, so it can play the role of $\mathbf{x}$ in your question.

Now, in practice most wireless channels fluctuate over time. This means that the vector $\mathbf{h}$ is not fixed, but changes randomly in time. If the symbol rate is not too high, then the channel will remain essentially fixed during the transmission of several symbols. This time is called the "channel coherence time". When the symbol time is much smaller than the coherence time, we say the the channel exhibits "slow fading".

Over a slow fading channel, the training sequence is repeated approximately once every coherence time. During one coherence time interval, the equalizer remains fixed.


Proper HW circuit to decode PWM signal using Arduino Mini

Discretize Amplitude or Time – that is the (signal processing question)

Seeing that you're using a microcontroller running at at leas 8 MHz (probably 16 MHz): It's way easier to measure the time a signal was high than to first average out the PWM  cycle in analog domain and then convert that analog voltage to digital.

If this was a sigma-delta ADC, which are pretty popular for most applications these days: Technologically, the conversion from analog voltages to digital numbers happens through counting the time a signal was high; I think the Atmega328p uses a successive approximation ADC, though.

The technological reality is that for signals below a couple of GHz, it's easier to build fast digital silicon logic than to build exact analog silicon. Thus, the less complex, more accurate and often less power-hungry solution is always to rather measure in time direction than in amplitude direction.

Implementation

So, I'd stay with the 220 Ω resistor: It really limits the current in the misconfiguration case, and is negligible compared to the proper input impedance.

Whether or not to add a capacitor really depends on what amount of noise you're expecting: you'll need a lot of noise to convert a stable TTL "high" to a "low", or vice versa. It still makes sense to do a bit of filtering, but not really with 0.1 µF:

----|220 Ω| ----·------ GPIO PIN
       R        |
                = C?
                |
                v GND

If we understand the RC combination as low-pass filter, then 220 Ω -- 100 nF leads to a cutoff frequency of ca 8 kHz, so not enough to average out a PWM cycle (if you wanted to go the analog voltage route), but maybe a little too much if you wanted to the low/high edges nice and sharp for detection. Assuming no significant inteferers between that frequency and 80 kHz, I'd maybe start with 10 nF – I think the ATmega328p has smitt trigger inputs, so you're really pretty safe from transients.

So, you only need to implement a way of measuring time: for the large time scales you're talking about, this could easily done by waiting for the pin to become high, and then just counting the CPU cycles until it turns low again. But that blocks the CPU for anything else – it'd be more elegant to use to use the built-in Timer/Counter units of the microcontroller, which have been invented for exactly that reason. You can program them to start counting at the rising edge of an external input signal, and stop on a negative edge.


A complex signal can be written as $s(t)=a(t)+jb(t)$, which indeed corresponds to two real signals ($a(t)$ and $b(t)$) "embedded" in the single signal $s(t)$. What makes this interesting is that $a(t)$ and $b(t)$ can be easily recovered from $s(t)$. This is true whether you write the complex signal in rectangular or polar coordinates.

Now, regarding how to interpret the phase of a signal. Let's say that the Fourier transform of $s(t)$ is $S(f)=|S(f)|e^{j\phi(f)}$. Imagine that you have an infinitely narrow band-pass filter with center frequency $f_0$. Then, you connect the signal $s(t)$ to the filter's input and see the output in an oscilloscope. You'll see a sinusoidal wave, with frequency $f_0$, amplitude $|S(f_0)|$ and phase $\phi(f_0)$. Here, "phase" indicates the difference, in radians, between a zero-phase cosine $\cos(2\pi f_0t)$ and the sinusoid you see in the scope.

In that sense, the phase $\phi(f_0)$ can be thought of as a "delay"; the sinusoid in the scope is just a cosine with a delay, $\cos(2\pi f_0t+\phi(f_0))$.

As an application, consider a signal that goes through an LTI system. The system doesn't touch $|S(f)|$ but it modifies the phase $\phi(t)$ -- that is, it changes the phase of each individual sinusoid. If this results in the same delay for each and every sinusoid, then the system output will be a delayed version of $s(t)$. If some of the sinusoids suffer different delays, then the signal at the system's output will be a distorted version of $s(t)$.


help with bandwidth diagrams for digital transmission

The $x$-axis in the diagrams is normalized frequency (normalized by the data rate). Note that they define $S$, which is the maximum baud rate, and $S_{ave}$, which is the average baud rate. You get the maximum baud rate if your data is such that your signal is changing at a maximum rate. For NRZ-L or Manchester that would be a sequence of alternating $0$s and $1$s. However, on average the baud rate is smaller than that maximum baud rate. This is what is meant by $S_{ave}$, which is less than $S=N/r$.

For baseband signals such as these, the (effective) bandwidth is the maximum frequency above which there are no relevant frequency components. Since the required bandwidth is proportional to the baud rate, and since both Manchester and differential Manchester codes have $r=\frac12$ whereas NRZ-L and NRZ-I have $r=1$, the Manchester codes need approximately twice as much bandwidth as an NRZ scheme.


In several application the mean is removed by using a high pass filter. A first order high pass filter has the following structure:
$$y[k] = x[k] - x[k-1] + \beta y[k-1]$$
or
$$y[k] = \beta x[k] - \beta x[k-1] + \beta y[k-1]$$
where $0<\beta<1$ determines the steepness of the high pass filter, larger $\beta$ leads to a steeper/sharper filter.

NOTE: Multiplying $x[k]$ and $x[k-1]$ by $\beta$ (as in the second equation) ensures that the gain per frequency is at most one. Due to the next normalization step, using this $\beta$ for $x$ is redundant.

Normalization of the signals energy can be performed with a recursive estimation of the signal energy:
$$ P_x[k] = \alpha y[k]y[k] + (1-\alpha) P_x[k-1] $$
where $y[n]$ is the zero mean signal and $0<\alpha<1$ is a forgetting factor. Small values for $\alpha$ lead to large time-constants of the recursive filter.

The normalized output output $z[k]$ is calculated as follows:
$$z[k] = \frac{y[k]}{\sqrt{P_x[k]}}$$

The signal $z[k]$ is a zero-mean signal with unit variance. You can scale your variance to a desired level $\sigma^2$ by applying the gain $\sqrt{\sigma^2}$ to signal $z[k]$.

NOTE1: A signal in the rang [0,1] always has a non-zero mean if the signal is not zero all the time.

NOTE2: In order to guarantee the signal levels to be in a specified range such as [-1,1], either a compression technique has to be applied or you have to know the original abolute maxima and minima such that you can scale to the desired range.

NOTE3: The variance in combination with the distribution only tells something about the majority of the samples.


How to segment with false edge in in image?

You could try "Non-maxima suppression" This is usually used on edge detectors to thin the edge to only one pixel. When there are lots of edges close together you basically pick the one with the largest gradient, and assume that is the true edge

Another technique that might be applicable is double thresholding. Some edge detectors (canny for instance) actually perform multiple detection's with different parameters.  "Strong edges" are found in the the higher threshold (and probably the lower one as well) weaker edges may only be found using the lower threshold, and are only "true" edges if they connect to one of the stronger edges. Maybe you could use your current edge detector as the weaker classifier, and then make one that is harder to satisfy. This might eliminate the false positives that are only detected by the weak detector (i hope that makes sense)

good luck


Why plotting the fft result of a signal from -f to f gives different result compared to 0 to f?

Can you advise:

What is the sample interval
What did you use to calculate the spectrum.

Edit 2:
The top one looks right, plotting from 0 Hz to 2*Nyquist. Not sure how
you calculated the spectrum, but this may be the convention
with that package.
If you prefer to present the spetrum from -ve Nyq to +ve Nyq you can't
just plot the same results with new axes.
You will need to recalculate
either with a function which uses that convention, or else do the DFT
yourself, or else you could just move the results above nyquist
and use these for the negative part.
Edit 3:
Nothing wrong, It is just an axis labelling thing.
The way your program has calculated the FFT, the answer is normaly interpreted
as 0Hz to 2*Nyquist.
Or if you prefer you could relabel it from left to right as:
(0-> nyquist), then (-Nyquist > 0)
If you want it to present as (-Nyq > Nyq) you need to rearrange the
answer. ("Roll" the right half of the array across to the left.)