Shift estimation using fourier transform

Assume the image is 256x256, The shift in spatial domain corresponds to a linear phase in the frequency domain.

u=-128:127;
v=-128:127;
u=repmat(u,size(u,2),1);
v=repmat(v',1,size(v,2));

linear_phase= exp(-2*pi*1i.*(((v).*dy)+((u).*dx)));

dx shift in $x$ and dy shift in $y$, and A equals linear_phaseso i can get a lot of equations to solve for dx and dy

If dx or dy is $0.5/128$ the shift is one pixel in $x$ or $y$


What is the Technique to Find Variance of Estimation Error

I don't understand the subscript $n$ notation, however, in the least squares problem that is given by:
\begin{equation}
{\bf{y}}={\bf{H}}{\theta}+\bf{n},
\end{equation}

where ${\bf{n}}\sim\mathcal{N}(\bf{0}, \sigma^2I_N)$ is a zero mean additive white Gaussian noise and $I_N$ is the $N \times N$ identity matrix, the maximum likelihood and the least squares estimators are equivalent and are efficient estimators (i.e. unbiased and attains the Cramer Rao bound) and are given by
\begin{equation}
\widehat \theta_{ML}=\left ( \bf{H^TH} \right )^{-1} \bf{H^Ty} .
\end{equation}

Now, let's calculate the first two moments of the estimator:

\begin{equation}
E\left \{ \widehat \theta_{ML} \right \} = E\left \{ \left ( \bf{H^TH} \right )^{-1} \bf{H^Ty} \right \} = \\
\left ( \bf{H^TH} \right )^{-1} \bf{H^T} \underset{=\bf{H}\theta}{\underbrace{E\left \{  \bf{y}\right \}}} = \theta,
\end{equation}
where $\theta$ is the true parameter value.
Thus, the estimator is unbiased and its Mean Square Error(MSE) will equal its covariance.

Now for the covariance:

Remember that if $\bf{z}=\bf{Ax}$ then $cov\left ( \bf{z} \right ) = \bf{A} cov\left ( x \right ) \bf{A^T}$.
Since the mean doesn't change the covariance in our problem $\bf{y} \sim \mathcal{N}\left ( \bf{H\theta}, \sigma^2I_N \right )$, thus,

\begin{equation}
cov\left ( \widehat \theta_{ML} \right )=\left ( \bf{H^TH} \right )^{-1} \bf{H^T} cov\left ( \bf{n} \right ) \bf{H}\left ( \bf{H^TH} \right )^{-1} = \\
= \left ( \bf{H^TH} \right )^{-1} \bf{H^T} \sigma^2 \bf{I_N} \bf{H}\left ( \bf{H^TH} \right )^{-1} = \\
= \sigma^2 \left ( \bf{H^TH} \right )^{-1} \bf{H^T}  \bf{H}\left ( \bf{H^TH} \right )^{-1} = \sigma^2 \left ( \bf{H^TH} \right )^{-1}.
\end{equation}

Note that $\bf{H^TH}$ is symmetric so it equals its transpose. Same is correct for its inverse.

I believe that this answers your question completely.


Frequency extracted from MATLAB's $\tt fft$ is not very accurate

MATLAB uses one based indexing, DFTs are zero based, so you have to shift your index value by one.

The bin index corresponds to the the frequency in units of cycles per frame.  Your freq is in units of cycles per second.  The sample count is in units of samples per frame and the sampling rate is units of samples per second. The bin index calculation should be:

\begin{align}
\mathrm{bin} &= \mathrm{freq} \cdot \frac{N}{f_s}\\
15 &= 10 \cdot \frac{1500}{1000}
\end{align}

Where the units are:

$$\rm \frac{cycles}{frame} = \frac{cycles}{second} \cdot \frac{\displaystyle\frac{samples}{frame}}{\displaystyle\frac{samples}{second}} $$

Bin 15 is in index 16 of MATLAB.  Since this is a whole number of cycles per frame, all your other bin values should be zero.

Hope this helps.

Ced

Follow up

Your calculation should have been:

$$\rm \frac{Frequency\ Bin \times Sampling\ Frequency}{Number\ of\ Samples} = \frac{15 \times 1000}{1500} = 10$$

Bin $-15$, aka bin $1485$, will also be non-zero and the complex conjugate of bin $15$.  Each should have a magnitude of $2.5$.  Your DC bin, aka bin $0$, should have a value of $-6$.  The corresponding MATLAB indexes will be $16$, $1486$, and $1$.


temporal high-pass filtering to approximate phase derivative

I am not familiar with the algorithm you have proposed, but I believe I can answer the question you asked.

Once you have computed the DFT of your image, what you get is a 2D array of complex number. Each value in this array represent the amplitude and phase of a specific spatial-frequency (not time-frequency). Basically, what your paper tells you, is that you must loop every 1024 value into your 32x32 array.

Consider the output of the 2D DFT as a new image. You need to loop every pixel, extract the angle of the complex number and filter that angle over time, and you do this 1024 times.

I believe your algorithm applies to a video feed, if that is correct, each pixel instensity is located in a 3 dimensions reference system : $I(x,y,t)$. After the 2D DFT, you get new values (complex values) that are located in a new 3D reference system : $I'(\omega_x, \omega_y, t)$. If you look at a single point (or a "frequency pixel", if you allow me the term)  in your 2D arrray individually, you have a time-series of complex values. You could denote it $I'_{pixel}(t)$. Consider the angle of each single point like this $\theta(t) = angle(I'_{pixel}(t))$. This is the time series that you must filter over time.

At this point any discrete filter could do, you'll have to figure which one gives you best result. The filter will be something like :

$$ \theta_{filtered}(t) = b_1\theta(t) + b_2\theta(t-1)+...+a_1\theta_{filtered}(t-1)+...$$

Where you choose the $a_n$ and $b_n$ coefficients to get a performant high pass filter.

And the nice part about matlab is that you don't have to iterate over each point as it is designed to do matrix operation. You could go with something like the code below, assuming your software is running in real time and that you grab a new frame everytime you call the function

img = grab_new_frame()
img_dft = fftshift(fft(img)) % 2D spatial DFT
angle_map = angle(img_dft)
filtered_angle = angle_map *b1 + last_angle_map*b2 + last_filtered_angle*a1 % Where an and bn are scalars
last_angle_map  = angle_map
last_filtered_angle = filtered_angle

Good luck!


How to Calculate the Local Gradient of an Image in MATLAB

Ok,
First of all, pay attention it is a calculation per pixel using a sata from a blog.
Basically summing the Gradient norm over a block / windows.

To calculate what you submitted above do the following:

1.Define the Finite Differences filters to approximate the gradient. Moreover build the "Win Summation" kernel which sums all elements within a window.

mLocalSumFilter = ones(filterWinSize, filterWinSize);
vGradXFilter = [-0.5, 0, 0.5];
vGradYFilter = [-0.5; 0; 0.5];

2.Calculate the image gradients:

mInputImageGradX = imfilter(mInputImage, vGradXFilter, 'replicate', 'same', 'corr');
mInputImageGradY = imfilter(mInputImage, vGradYFilter, 'replicate', 'same', 'corr');

3.Calculate the norm of the gradient vector. Afterwards sum the norm of the gradient within the window.

mInputImageGradNorm = sqrt((mInputImageGradX .^ 2) + (mInputImageGradY .^ 2));
mInputImageGradNorm = imfilter(mInputImageGradNorm, mLocalSumFilter, 0, 'same', 'corr') ;

This is it.
Enjoy.


The first and the obvious algorithm for the solution of this problem is the EM algorithm. This algoritm is getting worse in the number of parameters. It is also too much dependent on the initial conditions.

The drawbacks of this algorithm can be circumvented by parallel computing. This means you have $N$ EM algorithms run simultaneously which have different initialization points. The final algoritm takes all parameters estimated by N processors and based on the wellness of the fit gives the final decicion.

Another method is to use the central moments until Mth order. When written in a parametric form each xth order moment gives a linear/non-linear combination of parameters which are to be estimated. From here one can construct a test. Of course the moments should also be estimated from the data.

Another method would be to employ a better search algorithn. That is called particle swarm optimization. It will run similar to EM. However it is a stochastic algorithm have different characteristics than EM such as to be able to deal with more number of parameters. It is a kind of team optimization algorithm.

Thats all in my mind at the moment. Good luck with the experiments.


How to detect Premature ventricular contractions (PVCs) in a ECG-signal

I'm an audio guy and not a biomed guy, but if the ECG is decently periodic except for an occasional artifact or anomaly like a PVC, you can do "pitch detection" (estimating the period or fundamental frequency) of the waveform and tune a comb filter to null out the periodic component.  Then the anomalies will be the only thing left after subtracting or nulling out the periodic component.

Then you can perhaps have a dictionary of waveshapes for possible anolamies and cross-correlate (sliding correlation) each of the waveshapes against this remaining waveform.  That operation is, essentially, what we call a "matched filter", which is used to recognize "grains" of signals.


How can a stable system have a RHP pole?

It can't.

  Now, I have that if I do the root locus I have: (image left out of quote)
  so I have a pole in the RHP.

This is where you misunderstand.  The word "locus" means "The set of all points whose coordinates satisfy a given equation or condition".  The traces in that root locus show all possible locations for the poles -- and there's a lot of possible locations for poles inside the stability region.

  But I have that with an appropriate value of the gain of the controller we have a stable system. But in this case the pole should be in the LHP, not in the RHP.

In this case the poles are in the LHP.  You would find the gain where all the poles are stable (my intuition tells me you want the real part of your complex poles is in the neighborhood of $-2.5\mathrm{s^{-1}} < \mathcal{Re}(s) < -2.0\mathrm{s^{-1}}$), and then your system will be stable.


Low recognition rates using GMM for image classification

Are you trying to model the distribution of the SIFT descriptors using GMM? Are you using the actual 128-element descriptors?  If so, then you are trying to model a distribution in a 128-dimensional space, which is not likely to work well.

You need to reduce the dimensionality. If you insist on using GMM, then the simplest thing to try is to use principal components analysis (PCA) to reduce the dimensionality of your descriptors.  I would not use more than 10 dimensions.  Here is a paper, that tried kernel density estimation on reduced SIFT descriptors among other things.

Alternatively, you could try other algorithms, such as Bag Of Features, or the Pyramid Match Kernel.


Where to learn about "analog prototype filters"?

"Analog prototype" filters are well-known analog filters that have specific desirable properties. They include (but aren't limited to):

Butterworth filters: maximally-flat passband response
Chebyshev filters: sharper rolloff than Butterworth, equiripple in either the passband or stopband
Elliptic filters: even sharper rolloff, equiripple in passband and stopband
Bessel filters: maximally-flat group delay

These prototypes can be used to design digital filters that have approximately the same characteristics, for instance by using the bilinear transform.


What's wrong with my implementation of the Matching Pursuit algorithm?

In MP it is very important which dictionary is used. For example, if the atoms are orthogonal The problem will converge faster.

Another interesting example, where two entries of the dictionary are the same, you can wither find both corresponding values in the sparse signal as 1 and zeros or as 0.5 and 0.5. You can find any number of representation for such an example. Which one will the algorithm converge to?

In general, the picking of the dictionary $D$ is a very important point in the MP algorithm. I recommend you take a look in the following link as well as researching a bit more yourself. The Fourier matrix, for example, is a common choice for the dictionary.

To answer your question, while I could not find a problem with your implementation of the algorithm, the data you are using might be problematic.


wavelet decomposition for time series signal

A signal is just a 1D image. So if you can make it work for images, why not for signals? Pun aside, discrete wavelets combine multiscale smoother and differential operators, so they are been used as trend and singularity detectors for a while. Wavelet coefficients can be turned into feature vectors, directly or by condensating them across subband. For instance, wavelet subband coefficients are often believed to follow generalized Laplace-Gauss model distributions: $$A_{\alpha,\beta}\exp{\left(-\frac{|x-m|}{\beta}\right)^\alpha}\,.$$

The following coefficient histogram is for instance obtained from a stationary wavelet decomposition on a seismic signal (source: A Primal-Dual Proximal Algorithm for Sparse
Template-Based Adaptive Filtering: Application to
Seismic Multiple Removal)

And you can form a vector of subband features with estimated $\alpha_j $ and  $\beta_j $ at each scale $j$ for classification. Using stationary wavelets adds some shift invariance (see Difference between “equivariant to translation” and “invariant to translation”). Then, wavelet coefficients can be used as time-frequency anchors as in the Shazam music matching technique (An Industrial-Strength Audio Search Algorithm), or the deep-learning mimicking technique from the scattering wavelet transform developed by S. Mallat and coauthors.

Why does this make sense with (stationary)  wavelets?

Wavelets may act as multiscale change detectors, with a potential derivative behavior of wavelet functions at different scales, if you think about Haar (a discrete derivative) or the Laplacian of a Gaussian function with the Mexican hat wavelet,
They however can capture trends through the scaling function, since mother wavelets have vanishing moments,
They somehow can separate structured data from some noises, with orthogonality or whitening properties,
With stationary wavelets, the related features become integer shift invariant.


How to calculate the bandwidth of a continuous signal?

I suppose the correct expressions you meant are:

a) $5\sin(2t)$

b) $\sin(2t) + \sin(3t)$

Bandwidth is defined as the difference between the highest and lowest frequencies of a given signal ou system.

With this in mind, signal a) has one single frequency of 2 rad/s and so its bandwidth is 2-2=0 rad/s.
Similarly, signal b) has 2 frequencies: 2 rad/s and 3 rad/s. So, bandwidth is 3-2=1 rad/s.


Why is "$1/\mathrm{Hz}$" unchanged while converting$\mathrm{V}^2/\mathrm{Hz}$ to $\mathrm{dB}/\mathrm{Hz}$?

Say you have a measurement in unit $U$ (V$^2$) per unit $u$ (Hz). I do claim (although not in public) that the logarithm is a kind of unitless transformation, see for instance What is the logarithm of a kilometer? Is it a dimensionless number? In other words, it is somehow $0$-homogeneous, while a power is $2$-homogeneous: multiplying by a constant yields a constant to some power.

First, the logarithm is  kind of homogeneous to a $0$-th power: it is the integral of $1/x$ (power $-1$). The derivative of an affine function is a constant. Derivatives decrease powers by one unit, integrals increase them of the same amount. The logarithm somehow stands with constants between $x^{-1}$ and $x^{1}$. Additionally, the coefficient of variation of $x^p$ when $p\to 0$ is about:

$$ \frac{x^p-x^0}{p-0}\approx \log x\,.$$

Second, a dB scale is often relative: to some signal measurement, to some "absolute" value. The $\log$ argument is better unitless: something in U divided by some unit $1$ U, called U$_0$. Thus, computing $\log(\mathrm{U}/\mathrm{U}_0)$ could be considered.

So  converting U/Hz into dB$_\text{U}$/Hz, to me, only affects the numerator.


Recall from Physics class that density is, mass per unit volume. As an example, the approximate density of water is 1 g/cm3, so 100 cubic centimetres of water weighs 100 grams.

Power Spectral Density (PSD) is simply average noise power per unit of bandwidth. Power is expressed in units of watts. In telephony measurements, power is commonly expressed in logarithmic terms with respect to a milliwatt, i.e., dBm. The unit (volume) of bandwidth is
1 Hz. So PSD may be expressed as dBm/Hz, pronounced, “dBm per Hz.”

PSD is to Total Power as Density is to Weight. All transmission systems have finite power and bandwidth. Let’s consider an ideal system with a flat PSD signature.

At each frequency between f1 and f2 Hz, Let PSD to be a flat -10 dBm/Hz, or 0.1 mW/Hz.

Then

PSD (dB) = PSD(dbm/Hz) * Bandwidth(Hz)

where bandwidth can be defined as f2 –f1.

If you want to calculate the total power there are more info available in the reference.

reference: PSD techniques by Peter Walsh, NCE


High performance DSP processor with better on-chip memory

That depends a bit on the constraints: are you out of code space or data space or both?

Anyway, the typical way of dealing with space problems is to use external memory in conjunction with a fast access internal cache. An alternative would be something like this http://www.ti.com/lsds/ti/processors/dsp/c6000_dsp/c66x/overview.page with a few gigs of DDR attached.

Obviously external memory is somewhat slower than internal memory, so the trick is to move data efficiently into internal memory, so that it's there when you need it. In many cases the processor will do a pretty decent job in managing this automatically, especially if it's mostly code (not data) and most MIPS consumption is in loops.

If that's not good enough, there are certainly strategies to manually optimize code for cache efficiency as well.


How to adjust reverse range for Hue in cvInRangeS?

You can generate two different threshold images, one with hue in the range from 0 to 10 and the other with hue from 170 to 180. To get an image with hue in range both from 0 to 10 and 170 to 180, just add the two threshold images together.

In C++ code this would look something like:

inRange(imgHSV, Scalar(0, sMin, vMin), Scalar(10, sMax, vMax), imgThresh);
inRange(imgHSV, Scalar(170, sMin, vMin), Scalar(180, sMax, vMax), imgThresh1);

Mat imgThreshCombined = imgThresh + imgThresh1;

You could probably do the same thing by finding the threshold image with hue in the range from 10 to 170 and then just subtracting it from a white image.


Using k Nearest Neighbour & Dynamic Time Warping for keyword spotting

If you want to use the k-NN in conjunction with DTW, then it's quite easy. Assuming that you have the speech corpora of keywords you should do the following at the training stage:

Each keyword should have more than one utterance recorded
Extract the MFCC's for all the recordings and store them - these are effectively your "model data" - training samples.

Classification consists of:

Extract the matrix of MFCC's for the test sample
Calculate the DTW of the test sample with each and every MFCC matrix that was extracted and the training stage.
Store the final distance for each comparison. Now you should have a vector of total DTW distances (of a length equal to the number of the training examples):

test vs. class1_sample1: 10.2
test vs. class1_sample2: 20.3
test vs. class1_sample3: 15.1
test vs. class2_sample1: 3.1
test vs. class2_sample2: 2.5
...

For each of the classes (keywords), sort the distances (from small to large)
Pick the k top distances for each class and average them. These are the per-class memberships.
Predict the keyword by choosing the class with smallest distance.


You're right. That's why your anti-aliasing filter needs to be an analog filter. You need to restrict the bandwidth of the analog filter before sampling.

However, you can do part of the anti-aliasing filtering in the digital domain by first using a higher sampling rate and a simple analog filter, and down-sampling after applying a digital anti-aliasing filter. This is explained in the answers to this question. But in any case you'll need an analog filter first.


Without information on  $Φ$, you can obtain almost anything, since $\lambda Φ$ could be a valid CS matrix as well. Generally, one imposes structure contraints, such as unit energy for their rows or columns.

This being said, compressive sensing does not compress data in a strict sense, and energy is a poor measure of compressibility. Entropy and norm ratios could be more interesting measures.


Spectrum of windowed version of original continuous signal

If you don't understand the difference between the Continuous Time Fourier Transform (CTFT), the Discrete Time Fourier Transform (DTFT) and the Discrete Fourier Transform (DFT), now would be a good time to read about them. The very short version is that the DTFT yields the (continuous-valued) spectrum of a sequence (i.e., a sampled signal). The DFT computation results in a sampled version of the DTFT. To apply the DFT requires a finite number of samples (i.e., a time-domain window) whereas such restrictions are not placed on the DTFT in general. On the other hand, the CTFT deals with continuous time signals. There is a lot more to all of this, and I recommend you read more.

It is true that a complex-valued signal $x(t) = \exp\left(j \omega_0 t\right)$ and the Dirac delta function form a CTFT pair. I can agree with you that the spectrum of this signal is discrete (nonzero for a finite number of frequencies, in this case a single frequency).

After applying a rectangular window in time to the complex-valued signal $x(t)$, the CTFT is a frequency-translated sinc function centered at $\omega_0$ with a lobe width inversely proportional to the size of the time domain window. This shows us that the result of truncating a time domain signal with infinite support in time and a discrete spectrum in frequency can lead to a new time domain signal with finite support and a continuous spectrum. So, the answer is yes, the interpretation given in the question is true.


A priori spectral efficiency of OFDM scheme

Taking a look to the source code of the liquid-dsp project, and to the ofdmflexframegen.c file in particular, I found a partial answer to my question.

The spectral efficiency is:
$$\beta = \frac{R_b}{W}~~~~~~~~~~~~~~~~~~~~~~~~~~~(1)$$
where $R_b$ is the bit rate and $W$ is the OFDM bandwidth.

The bit rate is:
$$R_b = \frac{B}{T}~~~~~~~~~~~~~~~~~~~~~~~~~~~(2)$$
where $B$ is the data bits (of the decoded payload) and $T$ is the total transmission time.

Then, putting (2) into (1):
$$
\beta = \frac{B}{TW}~~~~~~~~~~~~~~~~~~~~~~~~~~~(3)
$$
and the product $TW = N$ is the number of samples of the transmission.

As for $N$, we have:

$$
N = (M + c)(N_{preamble} + N_{header} + N_{payload})~~~~~~~~~~~~~~~~~~~~~~~~~~~(4)
$$
where $M$ the total number of subcarriers, $c$ the length of the cyclic prefix.
$N_{preamble}$, $N_{header}$ and $N_{payload}$ are the number of OFDM symbols of the preamble, header and payload, respectively. Each OFDM symbol is comprised of $M + c$ samples.

Finally, the spectral efficiency for such a system would be:
$$
\beta = \frac{B}{(M + c)(N_{preamble} + N_{header} + N_{payload})}~~~~~~~~~~~~~~~~~~~~~~~~~~~(5)
$$

However, practical measurements of the spectral efficiency of my communication system are a little higher of what this formula predicts. It is like if the measured bandwidth was actually smaller than the predicted one. The difference increases with the payload size. I would appreciate any extra idea in this regard.

EDIT 1: Taking a deeper look to my system, it is possible that the difference between the spectral efficiency predicted by Eq. 5 and what I am measuring is due to the presence of null subcarriers. Null subcarriers at both sides of the OFDM spectrum are not measured, but they are included in Eq. 5. I wonder if I may just substract the number of null subcarrier from the term $M+c$ in Eq. 5. Does that make any sense?


Intuitively, why is windowing function a low-pass filter?

Windowing is not a filter.

Windowing is a multiplication of two signals in time (the input samples with the window function: $x_w[n] = x[n] \times w[n]$ ).

What you get in the frequency domain is a (circular) convolution of the transforms of your signal and the window.

This convolution in the frequency domain can be seen as the spectrum being low-pass-filtered, giving a resulting spectrum with less detail.

[In the case of a filter, you get convolution in time and product in frequency.]


Very first values of y[n] in an IIR filter

Typically, you would initialize the filter state (which in your example includes $x[-1]$, $x[-2]$, $y[-1]$, and $y[-2]$) with zeros. This is equivalent to assuming that the filter is causal.


How to remove the noise without destroying the main edge?

Use bilateral filter or anisotropic diffusion first.

The effect of anisotropic diffusion is as the following:

. The MATLAB code can be found here.

Here is its effect on your image:

Finally, non-local means is a also a good way to get rid of the noise. You might also want to take a look into that. I warn you though, it is slow.


How to extract the people and objects from an image separately?

I think what you're looking for is called "category level object recognition". As far as I know, it's still an active research subject. Meaning you can't expect to find high-quality libraries that you can use out of the box without understanding the magic behind. And even state-of-the art research code only achieves around 90% correct classification.

A common approach is the so called "bag of (visual) words" model. Meaning: You take a lot of labeled training images containing various objects. You use a keypoint detector&descriptor like SIFT or SURF to find the "visual words", i.e. similar-looking patches in the training examples, and you train a machine learning algorithm (like a support vector machine) to recognize the category labels based on the presence or absence of these words in each training image.

Then you again use the same keypoint descriptor, e.g. using a sliding window-technique on your test image, find the closest visual words, feed that information to your learning algorithm and you get an object classification for the sliding window (if you're lucky, that is).

Sadly, if you want to do that with Mathematica, you'll have to reinvent the wheel a lot. The SURF keypoint detector/descriptor is included in the ImageKeypoints function, but I don't think you can make it use sliding windows. And there is no machine learning functionality to speak of.


matching the mouth shape of an audio file(lip sync)

The simplest solution I can think of is to use an envelope follower, the higher the amplitude, the more open the mouth is. This is pretty simple but will one give tell you how 'open' the mouth should be, not the shape. This should already be quite effective for a cartoon. To make it more realistic you should probably try to detect vowels (Oh and Ah for instance) by measuring the formants. http://en.wikipedia.org/wiki/Formant gives a list of average frequency's for some vowels. Of course these are only average, so you need to experiment.

I don't think you need to recognize exactly which words are spoken, this will be much more complex, and there are quite a lot of words in each language, and quite a lot of languages... :-)

If you add consonant detection as well, you probably have all the information you need...


How is WSNR (weighted signal to noise ratio) defined?

There might be several versions, yet the most probable is as follows. An SNR is an energy-dependent measure. In the time or spatial domain, it is agnostic to location: all samples have the same "weight" (uniform weighting). As an energetic measure, it can be computed in the Fourier domain as well (or any orthogonal transformation).

Then, from the sensory system, we know that in hearing, or vision, not all frequencies have the same perceptual weight: a human hearing a loud 35 kHz tone won't probably not notice it, while a bat or a dog will notice. Hence, the relative importance of a sound should be weighted by the perceptual sensitivity of the subject (bat, dog, human). This yields a different weighting across different frequencies. One shall not care about a 35 kHz tone for a human ear, and give it a little weight, compared to the [1-2] kHz band. Same happens for vision. Such weightings are used for instance in JPEG compression.

So instead of computing, in the frequency domain some:

$$-10 \log \sum_k \frac{\left(\hat{s_k}-\hat{n_k}\right)^2}{\hat{s_k}^2}$$
with $k$ an index of frequency bins, you use a predefined weighting for each $k$:
$$-10 \log \sum_k w_k\frac{\left(\hat{s_k}-\hat{n_k}\right)^2}{\hat{s_k}^2}\,,$$
where  $w_k$ give less weight to less important frequencies. In sound, you can find A-weighting curves (and B-, C-, D- and Z-).

In images, a Matlab code function ratio = wsnr(orig, dith, nfreq) gives references to ancient papers:

J. Mannos and D. Sakrison, "The effects of a visual fidelity criterion on the encoding of images", IEEE Trans. Inf. Theory, July 1974.
T. Mitsa and K. Varkur, "Evaluation of contrast sensitivity functions for the formulation of quality measures incorporated in halftoning algorithms", IEEE Int. Conference on Acoustics, Speech and Signal Processing, 1993


Problems computing the DFT of finite length sequence


ValueError: operands could not be broadcast together with shapes (194,259,3) (600,800,3)

It sounds like the error is indicating that the two images aren't the same size.  So you can't subtract them point for point.


Actually, the Z transform is not really a proper transform, just a re-interpretation of the sequence of samples as coefficients of a formal Laurent series.

In some cases the formal Laurent series converges, if it does, it does so on an annular region in the complex plane. For useful signals (stable, summable, exponentially decaying) this annulus contains the unit circle, and the evaluation of the Laurent series on the unit circle corresponds to the Fourier series.

The interesting point of connecting a signal sequence to a periodic function on the unit circle is the inverse transformation, that many useful sequences are sequences of Fourier coefficients. And of course that convolution of signals corresponds to point-wise multiplication of the functions.


Using Chebyshev polynomials in wavetables

Chebyshev polynomials are mainly used for waveshaping, where an input signal is distorted by applying a memory-less non-linear function to it. This shaping function is stored in a look-up table (wavetable). If you use Chebyshev polynomials, then the polynomial just defines the shaping function, and the wavetable is generated just like for any other shaping function.

If the output is given by $y=f(x)$, where $x$ is the input, then the wavetable contains values $y$ for a given range of input values $x$. The function $f(x)$ can be defined in any desired way, and one way to define it is to use Chebyshev polynomials, which for a pure sinusoidal input generate a given harmonic.


OpenCV - What type of a result does cv2.goodFeaturesToTrack() return?

The cv2.goodFeaturesToTrack() is mainly used for detection of corners. See THIS PAGE

If you observe the type of the value returned by  cv2.goodFeaturesToTrack(), it is numpy.ndarray, which is a multi-dimensional array according to THIS DOCUMENT

So the function cv2.goodFeaturesToTrack() returns all the locations of the corners present in the gray scale image.

This is what a typical output would look like, the corners are marked in blue:


Lightweight test for bimodal distributions?

An easy way, especially one that is used in the detection and SNR estimation of PSK signals in a noisy signal, is based on stochastic moments of the received amplitude:

The kurtosis can be used as measure of how little "centered" a symmetric distribution is.

You could calculate a sample kurtosis frequency bin-wise in your spectrogram; those bins with high curtosis are likely to contain more than just Gaussian noise.

In a slightly more elaborate version of that algorithm optimized for detecting 2-FSK, you could use that to identify "candidate bins", subtract their mean power to have zero mean variables, and then calculate the covariance of pairs of these (i.e. calculate a covariance matrix). You're looking for a strong negative covariance – whenever one bin has above-average power, the other has low power.


Extracting fundamental amplitude/phase from only half a period of a pseudo-sinusoïd

If you know the frequency (and it seems that you do), simply construct a basis set of vectors from

1) DC (all ones)

2) Half a Sine wave

3) Half a Cosine wave

These form an orthogonal basis just like a DFT.  In fact you can think of it as a "one half bin and DC bin" specialized DFT.

Discard the DC reading (unless you care) and use the Sine and Cosine coefficients just like a regular DFT bin from which you can calculate magnitude (which gives you amplitude with a multiplication) and the phase (using atan2).

This will give you a best fit answer.  No window function and no zero padding.

Small correction (sorry):  The basis is not orthogonal.  You will need to solve a system of equations.

$$ x[n] = A \cos( \omega n ) + B \sin( \omega n ) + C $$

As a vector equation:

$$ \vec x = A \vec c + B \vec s + C \vec u $$

With some dot products, you get the system of equations with $A,B,C$ as the unknowns to be solved for.

$$ \vec x \cdot \vec c = A \vec c\cdot \vec c  + B \vec s\cdot \vec c  + C \vec u \cdot \vec c  $$

$$ \vec x \cdot \vec s = A \vec c\cdot \vec s  + B \vec s\cdot \vec s  + C \vec u \cdot \vec s  $$

$$ \vec x \cdot \vec u = A \vec c\cdot \vec u  + B \vec s\cdot \vec u  + C \vec u \cdot \vec u  $$

Note:

$$ \vec c \cdot \vec u = \vec u \cdot \vec c = 0 $$

$$ \vec c \cdot \vec s = \vec s \cdot \vec c = 0 $$

$$ \vec u \cdot \vec u = N $$

Where $N$ is your sample count.  Thus the equations simplify to:

$$ \vec x \cdot \vec c = A \vec c\cdot \vec c  $$

$$ \vec x \cdot \vec s = B \vec s\cdot \vec s  + C \vec u \cdot \vec s  $$

$$ \vec x \cdot \vec u = B \vec s\cdot \vec u  + C N $$

In order to get the best results with noisy data it is best to build the basis vectors over the full range of data you have and then solve the full set of equations.  A little more work, but likely worth it.

Okay, for the OP and maybe Ben.  Here is a sample program for the best fit of a partial cycle.  In this case 0.6 of a cycle.

import numpy as np
import matplotlib.pyplot as plt

#==========================================================
def main():

#---- Set Parameters, Base Frequency and Phase

        f_s = 50.0
        p   = 60.0 * np.pi / 180.0

#---- Set Domain Time and Angle for partial cycle

        t = np.arange( 0.0, 1.0, 1.0 / f_s )
        a = 2.0 * np.pi * t * 0.6

#---- Construct the Signal

        x = 10.0 * np.cos(        a + p ) \
          +  0.2 * np.sin(  2.0 * a )     \
          +  0.2 * np.sin( 10.0 * a )

#---- Build Basis Vectors

        c = np.cos( a )
        s = np.sin( a )
        u = np.ones( a.shape )

#---- Do Basis Dot Products

        cc = c.dot( c )
        ss = s.dot( s )
        uu = u.dot( u )
        cs = c.dot( s )
        cu = c.dot( u )
        su = s.dot( u )

        print cs, ss, uu, cs, cu, su

#---- Do Signal Dot Products

        xc = x.dot( c )
        xs = x.dot( s )
        xu = x.dot( u )

        print xc, xs, xu

#---- Build Coefficient Matrix

        M = np.array( [[ cc, cs, cu ],  \
                       [ cs, ss, su ],  \
                       [ cu, su, uu ] ] )

        print M

#---- Build Values Matrix

        V = np.array( [xc, xs, xu] )

        print V

#---- Solve for the best fit "Bin Value"

        R = np.linalg.solve( M, V )

        print R

#---- Calculate the Best Fit parameters

        A = R[0]
        B = R[1]
        C = R[2]

        D   = np.sqrt( A*A + B*B )
        phi = np.arctan2( -B, A )

        print D, phi, p

#---- Construct the Best Fit signal

        y = D * np.cos( a + phi ) + C

#---- Show the Results

        plt.plot( t, x )
        plt.plot( t, y )

        plt.show()

#==========================================================
main()

And these are the results:

These are the values:

Underlying target:

    10.0
p =  1.0471975512

Found best fit:

D   = 10.1939671316
phi =  1.04122684844

Changing the parameters a little bit for a smaller portion (0.3) of the cycle and moving it up by 20.

        x = 10.0 * np.cos(        a + p ) \
          +  0.2 * np.sin(  2.0 * a )     \
          +  2.0 * np.sin( 10.0 * a )     \
          + 20.0

Results:

Notice this isn't as accurate because of the smaller cycle portion.

D   = 10.4409622651
phi =  0.943824736628

However, if the signal is a pure tone:

        x = 10.0 * np.cos( a + p ) + 20.0

The results are still exact on less than a third of the cycle:

D   = 10.0
phi =  1.0471975512


How to prove $\int_{\Omega} \sum_{i=1}^{N} f_i(x)dx$ is equivilant with $\sum_{i=1}^{N} \int_{\Omega} f_i(x)u_i(x)dx$

Starting from your $E_2$ and exchanging the finite sum with the integral we get this:

$$E_2=\sum_{i=1}^N \int_\Omega f_i(x) u_i(x) dx = \int_\Omega \sum_{i=1}^N f_i(x) u_i(x) dx$$

Next we can use the disjoint cover $\{\Omega_i\}$ of $\Omega$ to split up the integral:

$$E_2 = \sum_{k=1}^N \int_{\Omega_k} \sum_{i=1}^N f_i(x) u_i(x) dx$$

The $x$ in the integral are now from $\Omega_k$ and we can apply your identity for the sum over $f_i(x) u_i(x)$:

$$E_2 = \sum_{k=1}^N \int_{\Omega_k} f_k(x) dx$$

This is almost your $E_1$, but not quite. We cannot exchange the integral with the sum because the integral depends on the summation index. Extending the integration region to $\Omega$ is also not possible without reintroducing a selection function. So your $E_1$ and $E_2$ are not equal. The closest thing to $E_1$ is what I have derived above.


In this line:

  anf->s.theta = -anf->p.gamma * anf->s.x1 * anf->s.theta * anf->s.err;

you seem to be conflating $\theta$ and $\dot{\theta}$ (the derivative).

Equation 5 of the paper is:

$$
\dot{\theta} = -\gamma x_1 \theta \left[ y(t) - \sum_{l=1}^n \dot{x}_l \right]
$$

Note that the left-hand term is the derivative of $\theta$, not $\theta$ itself.


%matplotlib inline is a magic command for IPython that allows you to add plots to the browser interface. It's not for the the terminal.
Sounds like you are using the Python interpreter at the command line, which isn't going to be useful for you as a beginner. Install jupyter notebook and use it instead of the command line interpreter to follow along with your course.
You do need to import IPython into your script (like any other module) if you want to use their functions during your script. This is not the same as installing a dependency from pip where it is downloaded and saved to your python local environment.


Normally a buffer of N points of data in the frequency domain covers, much like a two's complement number, -N/2 to +N/2-1.  So when you say that it doesn't cover N/2, it probably doesn't matter.  If it does matter and there is a "place" for N/2, just do as Paul suggested and set it to 0.


For a linear time-invariant filter, the output of the filter would be coloured noise.

Again assuming linear time-invariant filter, the effect of the filter on the output can be further characterized in terms of its power spectrum density (PSD).
Specifically, for an input $X(f)$ with power spectrum density $S_x(f)$ to a linear time-invariant filter with frequency response $H(f)$, the power-spectrum density $S_y(f)$ of the output of the filter is given by:

$$
  S_y(f) = |H(f)|^2 S_x(f)
$$

So, for a white noise input with unit power spectrum density $S_x(f)=1$, the corresponding PSD of the output would be
$$
  S_y(f) = |H(f)|^2
$$


Error in using Kalman Filter for 2D Position Estimation in Python

The main reason why your Kalman filter is not working is because you are not converting lat and lon values to kms. In the code below, I defined a new function called lat_lon_posx_posy which converts lat and lon values to px and py values in mts. You will need to make the following changes to your code.

Include the following function

import utm

def lat_log_posx_posy(coords):

     px, py = [], []
     for i in range(len(coords.lat)):
         dx = utm.from_latlon(coords.lat[i], coords.lon[i])
         px.append(dx[0])
         py.append(dx[1])
     return px, py

Change demo_kalman_xy() to the following:

def demo_kalman_xy():

    px, py = lat_log_posx_posy(coords)
    plt.plot(px[::18], py[::18], 'ro')
    plt.show()

    x = np.array([px[0], py[0], 0.01, 0.01]).T
    P = np.array(np.eye(4))*1000 # initial uncertainty
    result = []
    R = 0.01**2
    for meas in zip(px, py):
        x, P = kalman_xy(x, P, meas, R)
        result.append((x[:2]).tolist())
    kalman_x, kalman_y = zip(*result)
    plt.plot(px[::18], py[::18], 'ro')
    plt.plot(kalman_x, kalman_y, 'g-')
    plt.show()

Final Plot:


Multisampling of an ADC channel in a PWM cycle in TI C2000

Assuming your PWM above is ePWM1.
1 - Configure ePWM2 to be synchronized to ePWM1.
2 - Adjust ePWM2 period to be N times smaller.
3 - Now trigger the ADC sampling to ePWM2. Assuming you sample on the bottom and top of ePWM2 triangle, you will sample twice per ePWM2 period and 2*N times per ePWM1 period.


Removing Gaussian Noise from a Signal to Get Minimum Value

If you have no prior data on the signal of interest there is nothing to do actually.
The more prior you have the better you can do.
For instance, if the only information you have is the Bandwidth of your signal, the best you can do is apply an Band Pass / Low Pass Filter.
If you know a sparse representation of your signal it will be great, as you can remove all the other coefficients (Generalization of the LPF / BPF).
In case you know something about how smooth it is, you can apply Total Variation like smoothing.
Remark
When you write Gaussian Noise I assume you mean Additive White Gaussian Noise (AWGN).


Using morphological operations for color images processing

Take a look to this paper from Santiago Velasco, it's about Conditional Toggle Mapping, and it's application to denoising.


How do real-time convolution plugins process audio so quickly

Real-time low-latency partitioned convolution reverb with a long impulse response works by dividing the impulse response into unequally sized partitions. The shortest partitions (blocks) are at the beginning of the impulse response, and the partition length grows towards the end of the impulse response:

Each partition length can be processed separately, giving for each a latency that is for fast Fourier transform (FFT) -based convolution perhaps a couple of times the partition length. Convolution outputs for all partition lengths are summed forming the final composite output. When partitioning the impulse response, the partition of a particular length must not appear before the time that equals the latency of convolution of that partition length. A time-domain finite impulse response (FIR) filter can be used for the beginning of the impulse response because that may be more efficient than doing the low-length convolution by FFT.

What is the optimal partitioning for a given impulse response length and the optimal split point between time-domain FIR and FFT convolution should be tested on the target platform. A challenging aspect of non-uniform partitioned convolution is scheduling of the work.

Non-uniform partitioned convolution was popularized by William G. Gardner, Efficient Convolution without Input-Output Delay, J. Audio Eng. Soc., Vol. 43, No. 3, 1995 March, available online at http://www.cs.ust.hk/mjg_lib/bibs/DPSu/DPSu.Files/Ga95.PDF.


What is the most lucid, intuitive explanation for the various FTs - CFT, DFT, DTFT and the Fourier Series?

I wrote this handout as a complement to Oppenheim and Willsky. Please take a look at Table 4.1 on page 14, reproduced below. (Click for larger image.) I wrote that table specifically to answer questions such as yours.

Note the similarities and differences among the four operations:

"Series": periodic in time, discrete in frequency
"Transform": aperiodic in time, continuous in frequency
"Continuous Time": continuous in time, aperiodic in frequency
"Discrete Time": discrete in time, periodic in frequency

I hope you find these notes helpful! Please feel free to distribute as you wish.


Precalculate the Normalization of a filter size

Each pixel of the filter has a magnitude (intensity). The square function equivalently calculate the power density of the pixel. The normalization is to scale the filter power to 1. That's where I^2 at the denominator comes from.


Measurement of moving objects using computer vision

how difficult is it to stitch multiple images based on speed to get a complete object?

Trivial. If you know the speed and the time at which an image was taken, it's just a linear translation...

  My first instinct is to use line detection algorithms for detecting object edges and calculate distance between them, based on speed and time.

Excellent approach!

Just: reduce your image to a single row, maybe the average of the 10 center rows from your picture.

That reduces your problem from 2D to 1D.

Then, the question you have is "when does that start and stop being black?", which is trivial to answer if you just convert that color row to its brightness, and use a threshold.


Is wavelet analysis useful for 1D signals?

They absolutely are applicable. They are used for signal compression, denoising, feature extraction, filter-bank signal processing and many more things. In fact, a course of wavelet analysis will (should) stress one-dimensional analysis, which 2D is an extension of.


Is the incoming sample in a standard DPCM already quantized?

No it is assumed that m[k] is not quantized, as the main purpose of DPCM is to quantize the difference between samples rather than quantize the whole samples.

ofcourse if m[k] is an image it must been quantized as the video sensor in any digital camera is a digital sensor by definition so it is doing some means of quantization.


Difference between $\mathbb{E}[\mathbf{x} \mathbf{x}^{\rm{H}}]$ and $\mathbb{E}[(\mathbf{x}-\boldsymbol{\mu}) (\mathbf{x}-\boldsymbol{\mu})^{\rm{H}}]$

For the first case, as you wrote, it means the elements are not correlated. Since this is a Gaussian Random Vector it means the elements are independent.
It means that at most only one element of $ \boldsymbol{\mu} $ is not zero.
Since if there were more than 1, the matrix $ \boldsymbol{R} $ wasn't diagonal.

Update

Let's define $ \hat{\boldsymbol{x}} = \boldsymbol{x} - \boldsymbol{\mu} $, namely its centered version.

Looking ta the correlation matrix elements:

$$\begin{align*}
\boldsymbol{R}_{i, j} & = \mathbb{E} \left[ \boldsymbol{x}_{i} \boldsymbol{x}_{j} \right] = \mathbb{E} \left[ \left( \hat{\boldsymbol{x}}_{i} + \boldsymbol{\mu}_{i} \right) \left( \hat{\boldsymbol{x}}_{j} + \boldsymbol{\mu}_{j} \right) \right] \\
& = \mathbb{E} \left[ \hat{\boldsymbol{x}}_{i} \hat{\boldsymbol{x}}_{j} \right] + \boldsymbol{\mu}_{j} \mathbb{E} \left[ \hat{\boldsymbol{x}}_{i} \right] + \boldsymbol{\mu}_{i} \mathbb{E} \left[ \hat{\boldsymbol{x}}_{j} \right] + \boldsymbol{\mu}_{i} \boldsymbol{\mu}_{j} \\
& = \Sigma_{i, j} + \boldsymbol{\mu}_{i} \boldsymbol{\mu}_{j}
\end{align*}$$

Now it is really easy to see what is needed for this element to vanish.


Processes/Transforms involved to get brainwave data from raw EEG? (Autocorrelation confusion)

And the power spectral density is the Fourier transform of the autocorrelation.
I don't know about EEG data, but you can compute autocorrelation either in time domain or frequency domain. For instance if you need only a few samples of the autocorrelation it may be computed more efficiently using inner products than it would by $FFT^{-1}(|FFT(x)|^2)$. You know that the autocorrelation and the power spectral density are a time-frequency pair, so if the nature of information harmonic responses, oscillation modes, (power spectral density). But if the nature of the signal is chaotic with some dispersion (memory effect) you use the time domain (autocorrelation). autocorrelation are very good to analyze channel response in communication, and the signals in the brain are communication signals, so to me it makes sense to use that instead of frequency domain that could be used for identifying vibration in an engine for instance.


Pitch perception - Why does our ear not directly detect the missing fundamental frequency?

Here is roughly how it works.

A typically musical note is made out of harmonics and the fundamental. The fundamental is also the spacing between the harmonics.
For 100 Hz fundamental you would get 100Hz, 200Hz, 300Hz, 400Hz, 500 Hz, etc.
For 200 Hz you would get 200Hz, 400Hz, 600 Hz, which is simlar to 100Hz but not the same. Note that specifically 300Hz and 500Hz are missing. The spacing here is 200 Hz
In practice you  rarely get all the harmonics, some are louder than others, some are missing all together, however, you still can easily hear the pitch
That's because you determine pitch by the spacing of the harmonics, not by the lowest.
For example, if you generate a sound by adding 500, 600, 700 &  800Hz, you will hear the pitch corresponding to 100Hz, although no frequencies below 500 Hz are in the signal at all. 100 Hz is the periodicity, not necessarily the frequency content.
The organ of corti does actually a mechanical Fourier transform. The excitation through the oval window creates a bending wave on the basilar membrane. Due to mechanical resonance of the membrane and some active feedback mechanism controlled through efferent nerves and hair cells, there is a strong localized peak in the bending wave, the position of which corresponds with the frequency
In our example above you would see peaks at 500, 600, 700 & 800 Hz positions, but NOT at 100 Hz. There is no 100 Hz in the signal, so there is no excitation at the 100 Hz location.
The actual pitch detection is done by identifying correlated spectral parts parts of the signal and determining if there is a constant spacing. The spacing is the pitch. I don't know of the top of my head where in the periphery that happens but I'm guessing it's most like down through some correlation mechanism


What mathematical function could be used to simulate this sound effect?

I guess you want to synthesize that sound, i.e. create a synthetic signal to be as close as possible to the original sound. Potentially, you can do this by creating a superposition of sine waves with varying amplitudes and phases. So you have to know the parameters for those sine functions.
From what I can see and hear, the sound has two signals, which are very similar (both visually and to ear). Moreover, if you look at the frequency spectrum, you will see, that they are quite similar as well. However, the spectrum seems quite populated with different frequencies, which will make it difficult to reproduce as a synthetic function. At least if using only the provided original sound file.
Here is the spectrum plot:

And here is a python code to make those observations. (The file with sound is downloaded to the folder of the python script.)
import matplotlib.pyplot as plt
import soundfile as sf
import sounddevice as sd

# get file
filename = "SoundEffect.wav"
x, f = sf.read(filename)
x1 = x[:,0]     # the first signal
x2 = x[:,1]     # the second signal

# play and show the original sound in time domain
sd.play(x, f)
status = sd.wait()
plt.plot(x)
plt.show()

# play the first signal
sd.play(x1, f)
status = sd.wait()
# play the second signal
sd.play(x2, f)
status = sd.wait()

# show the signals in frequency domain
plt.magnitude_spectrum(x1, Fs = f)
plt.magnitude_spectrum(x2, Fs = f)
plt.show()

Bear in mind, I am not familiar with how sound effects are done. Most probably, this sound effect is a more simple superposition of certain signal shapes. However, unless they are already known, it is hard to "reverse-engineer" them from the sampled version of the signal.


Precision measurement of sine wave amplitude with ADC

A sine wave has infinitesimally little bandwidth. By rotating, filtering appropriately and decimating, you can reduce the sample rate very much.

Each of these filtering operations is typically a summing operation, in which you "average" out noise (which isn't your main concern), but also get a more precise estimate for the amplitude.

Decimation in DSP is very commonly done. You'd probably want to do that anyway – 320 MS/s is really no fun to deal with, and you don't need that bandwidth.

Of course, you can also correlate with a synthesized sine, and measure the correlation coefficient to get the power/amplitude.

Other options are things like proper spectral estimators – there's a lot to choose from, including Welch's method, or Pisarenko-based approaches. Especially if your signal is noisy, these might be interesting, but it really depends a lot on what exactly you're measuring and how you model your noise.


The precise meaning of 'coefficient' in signal processing

There's really not much I can clarify for you here:

A coefficient really is just a factor in series. Whether that series is e.g. a Fourier series, a sum over powers of 2, filter delay taps or just a series of sampling instants can't be said. Context!


Is it right to talk about dynamic range on frequency domain?

the dynamic range is the difference between the largest and the smalles value of a quantity can asume, so yes, you can have a dynamic range for the magnitud of your frequency components, I'm not aware if there's another terminology when dealing with the frequency domain. One thing to consider is not only the dynamic range of your signals but also the dynamic range of the variables you'rr using to store it and to calculate the frequency response, altough in the frequency domain (and specially when working with dB) it's not a big issue.


Calculate the magnitude and phase of a signal at a particular frequency in python

You can find the index of the desired (or the closest one) frequency in the array of resulting frequency bins using np.fft.fftfreq function, then use np.abs and np.angle functions to get the magnitude and phase.
Here is an example using fft.fft function from numpy library for a synthetic signal.
import numpy as np
import matplotlib.pyplot as plt

# Number of sample points
N = 1000

# Sample spacing
T = 1.0 / 800.0     # f = 800 Hz

# Create a signal
x = np.linspace(0.0, N*T, N)
t0 = np.pi/6   # non-zero phase of the second sine
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(200.0 * 2.0*np.pi*x + t0)
yf = np.fft.fft(y) # to normalize use norm='ortho' as an additional argument

# Where is a 200 Hz frequency in the results?
freq = np.fft.fftfreq(x.size, d=T)
index, = np.where(np.isclose(freq, 200, atol=1/(T*N)))

# Get magnitude and phase
magnitude = np.abs(yf[index[0]])
phase = np.angle(yf[index[0]])
print("Magnitude:", magnitude, ", phase:", phase)

# Plot a spectrum
plt.plot(freq[0:N//2], 2/N*np.abs(yf[0:N//2]), label='amplitude spectrum')   # in a conventional form
plt.plot(freq[0:N//2], np.angle(yf[0:N//2]), label='phase spectrum')
plt.legend()
plt.grid()
plt.show()

And here is a useful manual with detailed explanations: reference.


Will the phase of the reconstructed data using a single MRI k-space point show any curvature?

Some courses tell that each point in k-space (spatial-frequency domain) represents a certain stripe pattern in spatial domain.

It does. A simpler way to think about what sort of "stripe pattern" would be generated by a point $(m,n)$ in k-space is to imagine a wave that propagates along the direction of $(0,0) (m,n)$. That wave has a frequency that is proportional to the distance of $(m,n)$ from the centre $(0,0)$. Also, you might find this and this useful too.

Consequently, if I reconstructed the image using only a single point in k-space data, the phase of the constructed image will change periodically. This means there will be only straight lines in phase map.

This is not strictly accurate, or maybe it is described in a vague way here. The "phase of the constructed image will" not "change periodically". It will vary linearly.
If the phase varied periodically, the reconstructed image would look like frequency modulation.
Phase is the first derivative of frequency.
Here is an illustrative example in Python (but easily generalisable in other platforms):
import numpy

Fs=100 # Sampling frequency in Hz
T=1 # Duration of signal in seconds

t = numpy.linspace(0,T, Fs*T) # Time vector (in seconds)

p = 2*numpy.pi*t # Phase vector

Here, p is our phase vector. As you can see, it is just a line that increases at a specific rate. If you feed this phase to a sinusoid, you get a specific frequency.
If you want to generate a higher frequency you make the slope that the phase varies steeper. For example y = numpy.sin(10*p). This will give you a 10 Hz sinusoid running for 1 second at a sampling frequency of 100Hz.
By saying that "...the phase [...] will change periodically..." what you imply is that p somehow oscillates. But what will happen if we get an oscillating p? (Frequency Modulation).
Here is a very simple example (it comes as a continuation of the above one):
from matplotlib import pyplot as plt

plt.plot(numpy.sin(20*p+6*numpy.sin(p)));
plt.xlabel("Discrete time")
plt.ylabel("Amplitude")
plt.show()

Which will give you something like:

You can see here that we start at a high frequency, we go to a lower frequency and then back up again to a high frequency. The phase here still "evolves" at a linear rate but on that ramp, we have superimposed a sinusoid to show what happens when "...phase change[s] [...] periodically...".
What you are asking, is actually answered in the Stanford video you cite within the first 5-10 minutes.
To understand the sort of distortion that the gradient brings to the image, you have to understand how does frequency map to space in the MRI.
The key to understanding that is the Larmor frequency. The Larmor frequency is what allows us to map frequency to space (or more accurately space to frequency) by varying the strength of the magnetic field.
In other words, by having a linear magnetic field, we map equal distances in physical space to equal distances in the k-space.
If the magnetic field varies non-linearly (first 5-8 minutes of the Stanford video), then where you excite (or receive from) is not where you think it is.
Therefore, equal distances in physical space do not map anymore to equal distances in the k-space.
This is what distorts the image.
An extreme example in 1 dimension is (again, as a continuation of the above)
plt.plot(numpy.sin(p-numpy.sin(0.25*p)));
plt.xlabel("Discrete time")
plt.ylabel("Amplitude")
plt.show()

This will give you something like:

This was supposed to be a "normal" sinusoid cycle but due to the fact that the phase does not change linearly, what we see here is a distorted version.
Here is the difference between the linear evolution of phase and the profile that causes this distortion:
plt.plot(p);
plt.plot(p-numpy.sin(0.25*p));
plt.legend(["Linear", "Non-Linear"])
plt.xlabel("Discrete phase")
plt.ylabel("Amplitude")
plt.show();

Which gives:

Again, the key to understanding what is going on here is to relate physical space to frequency and field strength through the Larmor frequency. If the gradient is not linear then equal parcels of physical space are not mapped to equal parcels of k-space and that is where distortion is coming from.
Hope this helps.
EDIT:
To the points raised:

...my question is whether the phase is a straight line or curved?

You seem to be referring to the phase spectrum as obtained from the Discrete Fourier Transform (DFT).
This has nothing to do with the reason behind the image distortion and impossible to "probe" from a single image.
The distortion of phase because of gradient non-linearities is something that happens BEFORE the data is captured. By the time the signal is captured it is distorted already.
Therefore, if you do not have a reference image (that is, an image that is distortion free), it is impossible to "see" the phase distortion. You could try infering it too, but that would be valid up to a point and for assumed distortion profiles. The only reason why we are able to tell that an image is distorted is because we know the anatomy.

I just want to reproduce the results in this video.

The results presented in the video come from a simulator. The very first equation presented at the top left is the one that is used to map gradient phase to k-space position. If you don't have $\vec{G}$, it is impossible to know how $k$ gets deformed for a given image.
You can remap the phase yourself with a function (that is, change the balance between the imaginary and real components of an image's complex DFT) and proceed with the reconstruction, or, you can achieve the end result of the distortion via an image geometric transformation. This transformation is literally a remapping of the pixel positions to new positions which is exactly what happens when the gradient contains non-linearities. For more information you might want to see this.


Is the sum of all FFT bins from Audacity equal to RMS amplitude?

This is expected and a consquence of Parseval's theorem.

Loosely speaking, it's a flavor of energy conservation: the total energy of the signal doesn't change when you transform it from the time domain into the frequency domain so the energy calculated in either domain must be the same.


Guitar Tuner or Guitar effects Algorithms?

Best book on the topic is probably http://www.amazon.com/DAFX-Digital-Udo-ouml-lzer/dp/0470665998.

There is some freeware stuff (very mixed quality) here: http://www.musicdsp.org/

Many companies spend a lot of effort in making good sounding effects and well working tuners to make money with it, so access to "good" non-commercial sources if fairly limited.


Image interpolation is not (exactly) the same as upsampling. Suppose that you have a grid $\Lambda$ which is a $640 \times 480$ grid of pixels. Each pixel has a location $(x,y) \in \Lambda$. Each pixel also has an intensity $I(x,y)$. This could be a gray-scale value or an RGB triplet. Now suppose that you want to know the intensity of a location that is not centered on one of the pixels. For example, suppose that you wanted the intensity at location $(17.35, 188.42)$. In this case, you would use the neighboring integer pixel intensities in an appropriate interpolation scheme (nearest neighbor, bilinear, bicubic, Lanczos etc.) to get the sub-pixel intensity. This is image interpolation.

Image upsampling is when you want to go from a grid $\Lambda$ that is $640 \times 480$ to a larger grid $\Gamma$ that is for example, $1280 \times 960$ pixels. In such a case, a simple algorithm is create the new grid filled with zeros. Then, if there is a pixel value in the original image at $(x,y)$, then take its intensity I(x,y) and fill it in the upsampled image at $I'(2x,2y)$. The $2$ is the $r$ (rate) referred in your formula. Then, we use an appropriate kernel (the $h$ in your formula) to filter the image and complete the upsampling. NOTE: This is a simple explanation. Please see Gozalez, Woods and Eddins, Digital Image Processing for a better treatment.


frequency domain to time domain with magnitude and phase

Okay, I don't use MATLAB, but I think this Python example represents what you may be asking.

import numpy as np

A = 3 + 4j
k = 5
N = 128

Z = np.zeros(N, dtype='complex')
Z[k]   = A
Z[N-k] = np.conj( A )

z = np.fft.ifft(Z)

print( z )

This will generate a real valued pure tone (all the imaginary values are zero to the limit of precision) signal with 5 cycles in the frame.  The "mirroring" means that for a real signal the DFT is a conjugate mirror image, i.e. $Z[k] = Z^*[N-k]$.  It is common convention to use X, but I like Z too.


What is the name of these simple filter algorithms?

out[i] = out[i-1] + (alpha*(in[i] - out[i-1]));

From Wikipedia's article on Low pass filters:

  That is, this discrete-time implementation of a simple RC low-pass filter is the exponentially weighted moving average

  […]

  The loop that calculates each of the n outputs can be refactored into the equivalent:

for i from 1 to n
     y[i] := y[i-1] + α * (x[i] - y[i-1])

So, yes, this is a "exponentially weighted moving average". You can call it that, a single-pole IIR low pass with a decay of $\alpha$, or just a single-tap recursive low pass.

I'm happy I could help you with this, it wasn't very hard!


Does the determination of the ergodicity of a signal force any changes in methodology?

The concept of space as used in the statement ascribed to Wikipedia is presumably that of the ensemble of all possible signals, not space as in "Space: the final frontier...".  In the absence of parallel universes in which different signals from the ensemble can be observed simultaneously, it is rare for engineers to determine that a signal is ergodic by comparing time averages to ensemble averages.  More commonly, signals are assumed to be ergodic so that the values of time averages can be used as estimates of ensemble averages.  For example, an engineer takes $1000$ samples of a stochastic signal and plots a histogram of values.  This looks vaguely bell-shaped and so the signal is modeled as a Gaussian random process whose mean is the mean of the $1000$ samples etc.  Ultimately, the proof of the pudding is in the eating.  If this Gaussian model leads to predictions that differ wildly from observed results, the assumption of ergodicity might be discarded, or maybe $10,000$ samples, or $100,000$ samples might be tried (for example, predicting traffic at $5:30$ pm based on observations of traffic between $2:00$ am and $3:00$ am is rarely accurate, but
observations over several weeks might lead to a better model).


Here is a basic answer with some rapid experiments to show you what is possible to do with Morpholocical operations. Note that it is just a proof of concept, not a viable solution and you should improve it to better fit your problem.

First filter your image to detect edges. I used Canny but you may use a more appropriate one.

Then filter the BW result to keep only vertical and horizontal (hypothetical) lines by eroding the result with line structural elements. Here is the result for the erosion with a horizontal line of 5px:

Lastly sum each eroded images on its corresponding dimension. Here is the result of the sum of the previous image on the 2nd dimension:

From here you have to chose a threshold in order to infer the presence (or absence) of a collage. In the previous graph the sharp peaks indicate the position of horizontal boudaries of frames. A threshold value of 100 may be appropriate as it indicated that more that $1/3$ of the horizontal line is populated with edge pixels.


What is the difference of each frequency response of partial filters in a polyphase method

The short answer is the polyphase filter converts a low pass filter into a series of all pass filters each with a different time delay. So it is a series of delays at even fractions of the time between samples of the lower sampling rate of the polyphase filter. By getting an output of the same signal at different fractional delays, we can combine these to have a higher sampled version of the same signal. Polyphase means many phases, and the ideal implementation would be all-pass filters modifying the delay only- in good polyphase filter design we approach this ideal.

The longer answer is the polyphase filter is simply an array of filters formed by decimating the coefficients of a FIR filter, resulting in each filter output being identical to each other in their decimated (reduced) frequency response and each with a fractional delay of the other given the decimation relationship:

$$z^{-m}D\downarrow  = z^{-m/D}$$

Meaning a delay of $m$ samples at the higher rate will be a delay of $m/D$ samples at the lower rate after decimating by D.

This turns a low pass filter that would have a passband covering the frequency range after decimation into all pass filters, each with a different delay, hence the name poly-phase meaning many phases (it is actually many delays).

Consider the following 96 tap FIR filter with coefficients generated from MATLAB/Octave:

coeff = firls(96, [0 .2 .3 1], [1 1 0 0])

Which has this frequency response:

The time delay in samples is given by the slope of the phase vs frequency consistent with using $z = e^{j\omega}$ for $z^{-1}$: one cycle delay will have a phase that extends negatively to $\pi$ as the normalized angular frequency extends to $\pi$. Linear phase FIR filters (such as this one) have a delay that is half the number of taps, so in this case the filter delay is 48 samples, so the linear phase slope would be given by $z^{-48} = e^{-j48\omega}$ which is consistent with the phase slope we see in this filter's passband.

If we decimate this 96 tap filter into four separate 24 tap filters, the sampling rate for each filter will be at 1/4 of the rate above. Decimation is resampling and causes all the higher frequency spectrum to fold into the new sampling rate, while maintaining the original spectrum that was at the lower rate (if that is not clear, see this post The condition in order to not loose information after upsampling and downsampling). So if we decimated the impulse response of the filter (the coefficients), the frequency response below $0.25\pi$ rad/sample will remain the same other than an amplitude scaling by 4 or -12 dB, while the response above that will fold in, distorting the primary response if we didn't reject the higher bands sufficiently. For example, here is the frequency response if we take the first sample and every 4th sample after that in our original 24 tap filter, to create a 6 tap filter:

freqz(coeff(1:4:end));

So now the filter is nearly an "all-pass" filter modifying only the phase with a delay. It is passing all signals within it's new sampling frequency that is 1/4 of the original rate, and it's delay is 12 samples instead of 96: The phase at $\omega = \pi$ would be given by the phaes of $e^{-12j\omega}$; $12\pi  180/\pi = 2160°$, consistent with the maximum phase in the plot below.

If we select the next sample, and every 4 samples after that, we will end up with the same frequency response, since we are decimating the same filter, except the phase will be advanced by one sample at the higher sampling rate, or as we first explained $1/4$ sample at the lower output rate. Thus we would expect the phase slope to decrease to $11.75$ samples (with the associated smaller phase slope).

Superimposing all four filters given by coeff(1:4:end), coeff(2:4:end), coeff(3:4:end), and coeff(4:4:end) reveals that at least out to $0.8\pi$ radians/sample, these are now all pass structures with incrementally varying delay. The magnitude is completely flat out to $\omega = 0.8\pi$ which would be the usable range of our filter, we assume there is no signal of interest above $\omega = 0.8\pi$:

And the group delay (which is the negative derivative of phase with respect to frequency, and the time delay for linear phase filers) result is 12, 11.75, 11.5 and 11.25 samples:

The four delays are very conveniently at 1/4 samples of the higher 4x sampling rate! So each filter passes our signal within the passband with no distortion but provides 4 different outputs each at a time offset of the underlying waveform of 1/4 of a sample.

Given a signal exemplified by the sinusoid in the plot below, each color group would represent the output of one of the four polyphase filters.


Some toolkits like Kaldi can handle longer inputs. Others indeed require 10 seconds or so. It depends on the tools you gonna use.


The first definition of the power spectrum is the correct one (at least for causal signals, otherwise you need symmetrical integration limits). The second one assumes that the Fourier transform of $x(t)$ exists, which is generally not the case for random signals. The fact that the power spectrum and the autocorrelation function form a Fourier transform pair is no contradiction with the first definition of the power spectrum. This can be shown as follows:

$$S_x(\omega)=\lim_{T\rightarrow\infty}E\left\{ \frac{1}{T}\left| \int_{-T/2}^{T/2}x(t)e^{-j\omega t}dt \right|^2 \right\}=\\=
\lim_{T\rightarrow\infty}E\left\{ \frac{1}{T} \int_{-T/2}^{T/2}x^*(t)e^{j\omega t}dt \int_{-T/2}^{T/2}x(t')e^{-j\omega t'}dt'  \right\}$$

The inverse Fourier transform of the power spectrum $S_x(\omega)$ should equal the autocorrelation function $R_x(\tau)$:

$$\mathcal{F}^{-1}\{S_x(\omega)\}=\frac{1}{2\pi}\int_{-\infty}^{\infty}S_x(\omega)e^{j\omega \tau}d\omega=\\=
\lim_{T\rightarrow\infty}E\left\{ \frac{1}{2\pi T} \int_{-\infty}^{\infty}\int_{-T/2}^{T/2}x^*(t)e^{j\omega t}dt \int_{-T/2}^{T/2}x(t')e^{-j\omega t'}dt'e^{j\omega \tau}d\omega  \right\}=\\=
\lim_{T\rightarrow\infty}E\left\{ \frac{1}{T} \int_{-T/2}^{T/2}x^*(t) \int_{-T/2}^{T/2}x(t')\left[\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{j\omega (t+ \tau-t')}d\omega\right]dtdt'  \right\}
$$
The expression in brackets equals $\delta(t+\tau-t')$, which gives

$$\mathcal{F}^{-1}\{S_x(\omega)\}=
\lim_{T\rightarrow\infty}E\left\{ \frac{1}{T} \int_{-T/2}^{T/2}x^*(t) x(t+\tau)dt  \right\}=\\
=\lim_{T\rightarrow\infty} \frac{1}{T} \int_{-T/2}^{T/2}E\{x^*(t) x(t+\tau)\}dt=\\
= \lim_{T\rightarrow\infty} \frac{1}{T} \int_{-T/2}^{T/2}R_x(\tau)dt=
R_x(\tau)\lim_{T\rightarrow\infty} \frac{1}{T} \int_{-T/2}^{T/2}dt=R_x(\tau)
$$

which shows that the chosen definition of the power spectrum is compatible with the fact that the power spectrum and the autocorrelation function are Fourier transform pairs. A similar derivation can be done for discrete-time signals.


Practical Least Squares FIR filter design: How to find a realistic desired response

In general you seem to be doing all the right things.

The linear phase term is required since an arbitrary magnitude and phase response is typically non-causal. That means it's impossible to fit with a causal filter and your fitting error is very high.

One way to estimate a "good" delay value is to sample your target on an FFT grid that's relatively dense and do an inverse FFT. Ideally this looks like an actual impulse response with a clear peak someplace and trailing off at the edge. Circulate this until you've captured the bulk of the energy.

The other options is to try brute force: just try a lot of delays and pick the one that gives you the smallest error.

Controlling the out of band behavior you can play around with the weights of the
"don't care" bands. It's probably best to reduce the weights until the out of band energy is just about the same as the in-band energy. If you want real bandpass behavior, just apply a regular bandpass on top of your least square design. Trying to get a decent high pass at 250 Hz will cost you a lot of coefficients otherwise.

Your target has a lot of detail at low frequencies which determines the minimum number of FIR coefficients. You could try to prefit the major low frequency wiggles with an IIR filter and then just use an FIR for the remaining difference to the target.


Can playing a wav file theoretically damage my speakers?

Generally yes, but only if the amplifier is oversized for the speaker or there are other design faults. Active speakers where the amplifier needs to be generously sized to deliver the expected momentary performance are usually protected by peak limiters and filter circuits.

For laptop speakers the probability of damage is extremely low.


Causal system, order of numerator and denominator

When $H(z)$ is rational, then the system is causal if and only if its ROC is the exterior of a circle outside the out-most pole, and the order of numerator is no greater than the order of the denominator.

A causal LTI system with a rational transfer function $H(z)$ is stable if and only if all poles of $H(z)$ are inside the unit circle of the z-plane, i.e., the magnitudes of all poles are smaller than 1.

so you also need to find the R.O.C.


How to train and test deep neural network using MFCC features?

You may want to try training your model on 2D arrays containing entire datasets (all 75 MFCC frames of each audio file).  You can try recurrent or convolutional early layers in your DNN to handle the 2 dimensional input shape.  You may need a far larger amount of training data.


FFT - radix 2 - decimation in time -frequency range

If you studied the materials I recommended you on comments to your previous question. You would notice that FFT algorithm takes advantage of certain symmetries to make computations faster. If you don't want to work with a whole frequency range simply select the appropriate sampling frequency and concentrate all your FFT resolution points there. If you can't choose your sampling frequency take the samples and then decimate  them. Alternatively you can also use the Goertzel algorithm to work only with the frequencies you select.


Does a 20 kHz sine wave sampled at 44.1 kHz need to be reconstituted before playback?

Reminder about Nyquist-Shannon

By the theory of sampling, sampling a signal in the time domain results in periodisation of this signal in the frequency domain, at the period Fs (where Fs is your sampling frequency). Shannon-Nyquist's theorem explains that in order to not loose information when sampling a signal, there shouldn't be overlapping of the periodic repetitions of the signal's spectrum. Considering that the period of repetition is Fs, the total width of your signal's spectrum should exceed Fs.

Moreover, analog signals are real which means their spectrum is symetric relatively to the f=0Hz axis. Thus, if your signal's FT must have a maximum width of Fs, it means that your maximum frequency to avoid overlapping should be Fs/2.

Reconstruction filter

Now that your signal is correctly sampled i.e. there is no overlapping of the spectrum patterns in the frequency domain, you want to reconstruct it in the analog domain for example for audio playback. As explained before, sampling results in a periodisation of the analog signal's spectrum. That means that to go back to a continuous signal, you have to extract only one of the periods of your signal. To do that without distorting your signal, you need an ideal brick-wall filter (infinitely sharp transition band) called a reconstruction filter.
Unfortunately, such a perfect filter is impossible to obtain. Indeed, its impulse response is a cardinal sine (Fourier transform of a gate) which by definition is infinite on both sides of the central lobe.

It can however be approached with very high order FIR filters or convolution with a windowed sinc (this is a pretty deep subject, there's plenty of resources on the internet).

Answer to the original question

Depending on the quality of your reconstruction filter, sampling your signal close to Nyquist will indeed result in a distorted signal. That doesn't mean Nyquist theory is wrong, it is only theory. Depending on your DAC, you often want to oversample your signal if you need very high fidelity. Again, sampling you signal at a higher frequency will result in a periodisation of the signal at a longer period, which basically gives more space for your filter's transition band. The wider the transition band, the more oversampling you will need. As an example, the linear interpolation that your plotting software uses has the frequency shape of a squared sinc. It will thus completely distort your signal.


Decimating a signal (selecting every Dth sample and discarding the rest) does not distort the signal within the passband in any way other than to cause aliases from higher frequencies to fold into the signal bandwidth.  Depending on how we model the system the phase may be effected since $z^{-n}$ is replaced with $z^{-n/D}$, but the phase will still be linear, so the signal itself will not be distorted. (In multi-rate control loop modelling applications for example I had to pay attention to this change in phase slope). For further information on this side detail see the last paragraph.

Therefore an "anti-alias" filter is required, no different than what is done in and A/D application, just in the digital domain (decimation is a digital-to-digital converter).

The ideal decimation filter would pass the passband with no distortion while completely rejecting all the energy in the image bands. If that could feasibly be done then the decimation operation itself would be perfect (no effect on the signal in the passband). Actual performance depends on how well we design the decimation filter.  So in your case, of a signal in the second Nyquist zone (that once sampled will be in the first Nyquist zone with spectrum inverted as you describe) will not go through any further distortion with a subsequent decimation (if linear phase filter is used the result will be linear phase).

Since the folding areas are in distinct bands, a multi-band filter is a good choice for decimation filter design (firls and firpm in Matlab support the design of multiband filters).

As far as the impact on noise due to decimation itself, see this post that further details how the variance of the noise is not effected by the decimation: Properties of Up- and Downsampling

However the devil is in the details, as with any digital design you must pay attention to quantization effects, filter realization and passband distortion etc but this is not due to decimation specifically.

To further explain what I touched on in the first paragraph regarding the change in phase slope, consider a unit delay prior to a decimate by 2. A unit delay (with transfer function $z^{-1}$) has a constant magnitude over all frequencies, and a phase shift that varies from 0 to $2\pi$ as the frequency goes from DC to the sampling rate. After the decimate by two, the transfer function of the unit delay becomes $z^{1/2}$ which has a constant magnitude over all frequencies with a phase shift that varies from 0 to $\pi$ as the frequency goes from DC to the new sampling rate. So in the digital domain where we work with normalized frequencies, the slope of the phase has changed, but with regards to the related analog frequency in Hz, the slope is the same: At our new sampling frequency in Hz we see the same phase shift we got at that particular frequency prior to decimation.


Finding out modulation index and DC offset

The sinusoid test signal is given by
$$
s(t) = V_0\cos(2\pi f_0)
$$
and its Fourier transform by
$$
S(f) = \frac{V_0}{2}\big[\delta(f-f_0) + \delta(f+f_0)\big].
$$
Inserted in $S_\mathrm{AM-DSB-C}$ it yields
$$
S_\mathrm{AM-DSB-C}(f) = \frac{AV_0}{4}\big[\delta(f-f_c-f_0) + \delta(f-f_c+f_0) + \delta(f+f_c-f_0) + \delta(f+f_c+f_0) \big] + \frac {Ac}{2} \big[\delta (f-f_c) + \delta (f+f_c)\big].
$$
Comparing that with the figure given we see that $AV_0/4 = 3$ and $Ac/2 = 4$. Additionaly, $f_0 = 10\text{ kHz, } f_c = 100\text{ kHz}$. Using $m = V_0/c$ and $V_0 = 6\text{ V}$ we find
$$
m = 1.5, c = \frac{V_0}{m}=4\text{ V}.
$$


Is there a relation between the dead zone of sigma delta modulator and filter order?

From the paper you link to:

  Modulators with a  higher order loop filter also exhibit a  dead-zone which can be calculated in a  similar manner. Because of the higher gain in the loop filter, the dead-zone of these modulators will be considerably smaller than the dead-zone of the first order SDM*. Very small signals (including system noise) will disturb the idle patterns of such SDMs, rendering the dead-zone undetectable in practical implementations.

And if you had a situation where system noise didn't disturb the idle patterns that cause dead zones, you could always inject some dither that's intentionally designed to have its spectral content at high frequencies.

* SDM = sigma delta modulator, in that paper's nomenclature.


Application of the time-shifting property in case of Fourier-Transform of cosine

Two things

Time shift equation is wrong. Should be $e^{-j \omega n_d}$.
You need to express the cosine as a time shift $cos(\omega (n - \pi /2 /\omega))$ Time shift is not the same as phase shift and the  time shift is actually pi/2 divided by the frequency


What is a low-pass filter for a sequence of tuples?

'Frequency' for a sequence of values is the measure of how sharp the values in that sequence change.

For example, the sequence $x_1(n)=[\begin{array}[cccccccc]&\cdots &1, &1, &1, &1, &\cdots\end{array}]$ has a frequency 0 but the  $x_2(n)=[\begin{array}[cccccccc] &\cdots &1 &1 &7 &1 &\cdots\end{array}]$ is a  higher frequency signal since there is a sharp change in sequence value (1 to 7 and 7 to 1).

A low pass filter is a system which attenuates the high frequency components. OR it reduces the sharpness (or introduces blur) in case of an image. OR a low pass filter does not allow sudden changes in values and smoothens these changes. Figure below shows low passing various sequences ($\style{color:blue}{blue:input}$ $\style{color:red}{red:output}$).

Consider a system whose output is a sequence $y(n)$ where n$^{th}$ value in $y(n)$ is the average of n$^{th}$, (n-1)$^{th}$ and (n+1)$^{th}$ value of input sequence is a simple low pass filter. If $x_2(n)$ is fed as input to this filter, then output, will  be $[\begin{array}[cccccccc] &\cdots &1 &3 &3 &3 &1 &\cdots\end{array}]$ and it is clear that the magnitude of change between consecutive values are reduced. or the high frequency components are attenuated.

Since you have a tuple, you have to separate it into two sequences ([x1,x2,....] and [y1,y2,....]). The low pass filtering should be done on these sequences separately.

In character recognition algorithms, low pass filtering is used in pre-processing to remove noises or unwanted small dots, having a width of one or two pixels. see an example of low pass filtering used for noise removal in image processing.


Why this two FIR filters have the same amplitude but different phase response

If you have a causal length $N$ FIR filter with impulse response $h[n]$, and with frequency response $H(e^{j\omega})$, and you invert it on the time axis, you get a new filter

$$g[n]=h[-n]\tag{1}$$

This filter is non-causal, but since it is an FIR filter, it can be made causal by shifting it to the left by $N-1$ samples:

$$f[n]=g[n-(N-1)]=h[N-1-n]\tag{2}$$

This is what you're doing. Now since

$$H(e^{j\omega})=\sum_{n=-\infty}^{\infty}h[n]e^{-jn\omega}\tag{3}$$

you have

$$G(e^{j\omega})=\sum_{n=-\infty}^{\infty}h[-n]e^{-jn\omega}=\sum_{n=-\infty}^{\infty}h[n]e^{jn\omega}=H(e^{-j\omega})$$

and

$$F(e^{j\omega})=G(e^{j\omega})e^{-j(N-1)\omega}=H(e^{-j\omega})e^{-j(N-1)\omega}$$

From (3) you can see that for real-valued $h[n]$ you have $H(e^{-j\omega})=H^*(e^{j\omega})$ (where $*$ denotes complex conjugation). So you finally get

$$F(e^{j\omega})=H^*(e^{j\omega})e^{-j(N-1)\omega}\tag{4}$$

and since $|e^{-j(N-1)\omega}|=1$, and $|H^*(e^{j\omega})|=|H(e^{j\omega})|$ you have

$$|F(e^{j\omega})|=|H(e^{j\omega})|\tag{5}$$

So the magnitude responses of the original filter and the time-reversed filter are identical. The phase responses are different: from (4) you have

$$|F(e^{j\omega})|e^{j\phi_F(\omega)}=|H(e^{j\omega})|e^{-j\phi_H(\omega)}e^{-j(N-1)\omega}\tag{6}$$

from which

$$\phi_F(\omega)=-\phi_H(\omega)-(N-1)\omega\tag{7}$$

follows. Since the group delay is the negative derivative of the phase, the group delays are related by

$$\tau_F(\omega)=-\tau_H(\omega)+N-1\tag{8}$$


I don't really understand the thing about cosines (as in: how is that helpful?) – a DFT is really just a mapping from a complex vector with $N$ elements to a complex vector with $N$ elements; and your calculation seems to be wrong, and I'm not sure where, but doing two of the three elements of the DFT manually might actually be enough to clear things up.

\begin{align}
X_0 &= \sum\limits_{n=0}^{N-1} x_n e^{-i2\pi\frac nN\cdot 0}\\
&= \sum\limits_{n=0}^{N-1} x_n e^0\\
&= \frac 14 (1+2+1)\\
&=1\\[2em]
X_1 &=\sum\limits_{n=0}^{N-1} x_n e^{-i2\pi\frac nN\cdot 1}\\
&=\sum\limits_{n=0}^{N-1} x_n e^{-i2\pi\frac n3}\\
&=\frac14 \left( 1 \cdot e^{-i2\pi\frac 03} + 2 e^{-i2\pi\frac 13} + 1 e^{-i2\pi\frac 23}\right)\\
&= \frac 14\left[1+ 2\cos\left(2\pi\frac13\right)+ i2\sin\left(2\pi\frac13\right)
+ \cos\left(2\pi\frac23\right)+ i\sin\left(2\pi\frac23\right)\right]\\
&\text{bit of basic trigonometry}\\
&=\frac14\left[
1+ 2\cos\left(2\pi\frac13\right)+ i2\sin\left(2\pi\frac13\right)
- \cos\left(2\pi\frac13\right) + i\sin\left(2\pi\frac13\right)
\right]\\
&=\frac14\left[
1+ 1\cos\left(2\pi\frac13\right)+ i3\sin\left(2\pi\frac13\right)
\right]\\
&=\frac14\left[
1-\frac12+ i3\frac{\sqrt3}2
\right]\\
&=\frac18+i\frac{3\sqrt3}8
\end{align}

Note how there's no taking of real parts or anything: the DFT of a sequence is complex in general, and you can't just drop the imaginary part. So, not quite sure why we're talking in control theory terms here, but no, you cannot just drop it. And: the discrete Fourier is not the same as the Laplace transform.


Having the impulse response allows one to simplify the input of a sampled linear time-invariant system down to one single unit-magnitude input sample.

It is usually fairly simple to build any time domain input from a single unit input (with a linear combination thereof), and thus, given the impulse response for that single input, the output for any time-domain-described input for that system.  Other basis vectors are possible, but may require more or less work (computational effort) for a given type of description of the input.


changing circular shift in hologram reconstruction

The problem is this line:

recon_image = ifftshift(ifft2(psi));

It should be

recon_image = ifft2(ifftshift(psi));

You need to shift the 0 frequency back to the corners before doing ifft2

I find fftshift / ifftshift and creating the frequency coordinates causes so many programming errors. This is my code for creating the coordinates (with 0 frequency in the corner):

function [ux,uy,r,th] = fftMesh(sizeV)
%FFTMESH - Create DFT cartesian and polar coordinates
%
% Inputs:
%
% sizeV     [width,height] of spectrum (same as image size)
%
% Outputs:
%
% ux        x in cartesian coords
% uy        y in cartesian coords
% r         radius in polar coords
% th        angle (rad CCW) in polar coords

rows = sizeV(1);
cols = sizeV(2);

[ux, uy] = meshgrid(...
    ([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
    ([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2))...
    );
ux = ifftshift(ux);   % Quadrant shift to put 0 frequency at the corners
uy = ifftshift(uy);

% Convert to polar coordinates
th = atan2(uy,ux);
r = sqrt(ux.^2 + uy.^2);

end

If you use this, you dont have to worry about using fftshift / ifftshift. Alternatively, just remove the ifftshift lines in the function to keep the 0 freq in the middle.

Updated code:

clear;

lambda = 632.8e-9/1.333; %wavelength
pix_size = 7.4e-6; %camera's pixel size
dx = pix_size;
dy = pix_size;

hologram = double(rgb2gray(imread('test.jpg')));
N = size(hologram);

tukey_win = tukeywin(N(1),0.075); %tukey window
% multiply image (single precision) with the tukey window to darken edges
trimmed_holo = single(hologram) .* (tukey_win * tukey_win');

prop_distance = 0.01:0.005:0.2; % propagation distances vector
fft_hologram = fft2(trimmed_holo); %Fourier transform of hologram

for k = 1:length(prop_distance) % do multiple propagation distances
% obtain transfer function
G = trans_func(lambda, N(1), dx, dy, prop_distance(k));

%multiply FT of impulse response times FT of hologram
psi = G .* fft_hologram;

% get slice of reconstructed optical field through inverse FT
recon_image = ifft2(psi);

%write image to file
imagesc(real(recon_image)); pause;
imagesc(imag(recon_image)); pause;
imagesc(abs(recon_image))

end

function G = trans_func(lambda, N, Delta_xi, Delta_eta, z)

M = N;
[m, n] = meshgrid(...
    ([1:M]-(fix(M/2)+1))/(M-mod(M,2)), ...
    ([1:N]-(fix(N/2)+1))/(N-mod(N,2))...
    );
m = ifftshift(m);   % Quadrant shift to put 0 frequency at the corners
n = ifftshift(n);

first_term  = (lambda^2 * (n + (N^2 * Delta_xi^2)  ./ (2 * z * lambda)).^2) ...
    ./ (N^2 * Delta_xi^2);
second_term = (lambda^2 * (m + (M^2 * Delta_eta^2) ./ (2 * z * lambda)).^2) ...
    ./ (M^2 * Delta_eta^2);
arg_root = 1 - first_term - second_term;
arg_root(arg_root<0) = 0;


Why there is error between My DFT and built in Matlab FFT?

Your error is the yellow line at the bottom.  i see no difference between your DFT and the built-in FFT from MATLAB.


Image Reconstruction:Phase vs. Magnitude

Figure 1.(c) shows the Test image reconstructed from MAGNITUDE spectrum only. We can say that the intensity values of LOW frequency pixels are comparatively more than HIGH frequency pixels.

Actually, this is not correct. The phase values determine the shift in the sinusoid components of the image. With zero phase, all the sinusoids are centred at the same location and you get a symmetric image whose structure has no real correlation with the original image at all. Being centred at the same location means that the sinusoids are a maximum at that location, and is why there is a big white patch in the middle of Figure 1.c.

The phase-only reconstruction preserve features because of the principle of phase congruency. At the location of edges and lines, most of the sinusoid components have the same phase. See http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/OWENS/LECT7/node2.html  This properly alone can be used to detect lines and edges, http://www.csse.uwa.edu.au/~pk/research/pkpapers/phasecorners.pdf , without regard to magnitude. So you can see that the phase information is most important.

Changing the magnitude of the various component sinusoids changes the shape of the feature. When you do a phase-only reconstruction, you set all the magnitudes to one, which changes the shape of the features, but not their location. In many images the low frequency components have a magnitude higher than the high frequency components, so phase-only reconstruction does look like a high-pass filter.

In short, phase contains the information about the locations of features.

You cannot add the phase-only and magnitude-only images to get the original. You can multiply them in the Fourier domain and transform back to get the original.


The word "Filter" in Image Processing word relates to "Neighborhood" based operation on pixels.

Filters in the context of Instagram and other Image Processing applications are combination of point operations (Works on single pixel, adjusting it values independent of its adjacent pixels) and filters which alters the pixel value according to its adjacent pixels.

If I were you I'd start with learning about Point Operations such as Brightness, Contrast, Hue, Saturation, etc...
Once you build some functions which alters those move on to filters.

Then you could combine them into "Effect".

For instance, Vignetting is easy to do:

Create a rounded mask. White out side a circle and black around it.
Smooth the mask using a Gaussian Blur.
Adjust the brightness of pixel according to the values of this mask.