How do the English pronunciations of "tiu" and "tou" differ?

First of all, don't expect an answer like "Compute such and such well-known function on the data and if the result is above such or such value it's tiuuu otherwise it's tou". You won't get anything robust or reliable with approaches like that.

The good representations of sound which can be used for this kind of discrimination task are either the Mel Frequency Cepstral Coefficients (MFCC) or the linear prediction coefficients (Compute the first 20 or so autocorrelation coefficients and use the Levinson-Durbin recursion to get the AR coefficients). Both are somehow independent of pitch (good), moderately robust to noise (good but can be better), but are still somewhat speaker and gender dependent. It's still a vast improvement over trying to guess formants from a spectrum and matching that to a per-phoneme table.

Does your application tolerate a training phase? If so, get your user to record a dozen of utterances of each word ; and use this as a training set for a simple bayesian classifier using a gaussian model. Use this classifier to discriminate new utterances.

Does your application need realtime processing? If not, you can get by with clustering. Record the whole "sentence", split it into "words" ; compute average feature vectors over each word ; and run K-means clustering to get 2 class out of it. One class will be "tou" and the other "tiu". The class whose "words" have the lowest average spectral centroid are likely to be the "tou".

A last word on speech recognition: if you have a vocabulary of two words, the cost of running a decent speech recognition toolchain will be that of:
* Computing the MFCC features (which you have to do anyway no matter what).
* Scoring a dozen of gaussians (which you have to do anyway).
* Finding the lowest cost paths in a FST with less than a dozen states.

This is quite pedestrian and nowhere near the complexity of natural language speech recognition. Since you have no experience in digital signal processing, doing so will probably be faster than whatever custom code you might write. Thus I very strongly suggest you to collect a small training corpus of utterances from people around you and start looking into "off the shelf" speech recognition toolboxes like HTK. Just using the built-in tools you will be able to train a model evaluate it on your data without writing code.


Plotting phase delay from first principles giving different result from phasedelay() function. Why?

The problem is the definition of the phase. The command angle() computes the phase $\phi(\omega)$ according to

$$H(e^{j\omega})=|H(e^{j\omega})|e^{j\phi(\omega)}\tag{1}$$

where $H(e^{j\omega})$ is the complex frequency response. At frequencies $\omega$ where the frequency response has zeros, the phase jumps by $\pi$. This is shown in the two left figures of the plot below. You can't compute the phase delay from $\phi(\omega)$ because of the jumps. Note that in your computation, the phase delay is correct until the first phase jump due to the first zero of the magnitude response.

The phase delay needs to be computed from a continuous phase $\varphi(\omega)$, which is implicitly defined by

$$H(e^{j\omega})=A(\omega)e^{j\varphi(\omega)}\tag{2}$$

where $A(\omega)$ is a real-valued but bi-polar smooth function. This continuous phase is shown in the top right figure below, and it is obviously linear. From this continuous phase, the phase delay can be computed by

$$\tau_{ph}=-\frac{\varphi(\omega)}{\omega}\tag{3}$$

and for the given linear-phase FIR filter the phase delay is constant:

$$\tau_{ph}=\frac{N-1}{2}\tag{4}$$

where $N$ is the filter length.

Note that the continuous phase $\varphi(\omega)$ cannot be obtained from $\phi(\omega)$ using the Matlab function unwrap(), because the latter just removes jumps by (multiples of) $2\pi$ due to phase ambiguity. Here we have to deal with actual phase jumps due to zeros of $H(e^{j\omega})$.


If I understand correctly, you want to verify the energy calculation in the frequency domain by computing the energy as

$$E_x=\int_{-\infty}^{\infty}|X(f)|^2df\tag{1}$$

with

$$X(f)=\mathcal{F}\big\{x(t)\big\}=\frac{A}{A+i2\pi f}\tag{2}$$

From $(2)$ we get

$$|X(f)|^2=\frac{A^2}{A^2+(2\pi f)^2}=\frac{1}{1+\left(\frac{2\pi f}{A}\right)^2}\tag{3}$$

With $(3)$ and with the substitution $x=2\pi f/A$, the integral $(1)$ becomes

$$E_x=\frac{A}{2\pi}\int_{-\infty}^{\infty}\frac{1}{1+x^2}dx=\frac{A}{2\pi}\arctan(x){\huge|}_{-\infty}^{\infty}=\frac{A}{2\pi}\cdot\pi=\frac{A}{2}\tag{4}$$


Signal Resampling by a non integer factor with convolution

You might be describing something similar to a canonical polyphase resampling filter.  If you can compute your filter+interpolation kernel on the fly (certain windowed Sinc’s and others) you can even do this for irrational sample rate ratios or factors.


Are SURF descriptors normalized with respect to orientation?

The SURF descriptor implicitly accounts for orientation of the interest points based on the distribution of edge intensities around it. So yes it is invariant to rotation. This is also mentioned in the description of the SURFPoints class:

  Orientation

  Describes the orientation of the detected feature. This value must be
  specified as an angle, in radians. The angle is measured from the
  X-axis with the origin at the point given by the Location property.
  Typically, this value gets set during the descriptor extraction
  process. The extractFeatures function modifies the default value of 0.
  Do not set this property manually. Rely instead on the call to
  extractFeatures to fill in this value. The Orientation is mainly
  useful for visualization purposes.

  Default: 0.0


I'll start with the Hough line transform.

The hough transform works by creating a buffer (2 dimensional in our case) which represents all the possible lines in the image.

Any possible line can be represented with those 2 numbers $(\rho, \theta)$

From wikipedia

Note this uses polar coordinates to avoid division by zero issues rather than your typical $y = mx + c$ for a line

The hough transform scans an image and casts votes into the 2 dimensional buffer corresponding to $(\rho, \theta)$ values.

The buffer slots that accumulate the most votes are likely to be attributed to the presence of a line for that $(\rho, \theta)$ value.

That's the basic idea of the Hough transform but it can get more complicated from here.

For selecting $(\rho, \theta)$ values from the buffer, methods can get more sophisticated, such as using kernel density estimation. Although I don't think opencv does that.

In the case of the Hough circle transform, we now have a 3 dimensional buffer.
$(x, y, r)$

Where

$x, y $:  coordinates of center of circle

$r $: radius

Now in practice there's more to making the Hough transform work than this.

Take a look at the source code for the icvHoughCirclesGradient function.

Its here

Also cvHoughCircles and cv::HoughCircles

you can see that the function is also doing canny edge detection and sobel filtering.

These operations have their own input parameters which you can also see as inputs to the hough function.

Examining the function signature

void cv::HoughCircles( InputArray _image, OutputArray _circles,
int method, double dp, double min_dist,
double param1, double param2,
int minRadius, int maxRadius )

Method: There seems to be only on method now, probably left this for future development

dp: this is the resolution of the votes buffer

min_dist : is the minimum distance between circles (ie: circles that are too close must be part of the same circle)

minRadius and maxRadius are just that.

param1 and param2 are a bit more cryptic

but if you spot the lines

int canny_threshold = cvRound(param1);
int acc_threshold = cvRound(param2);

You'll see these are related to the canny edge detection and the accumulator threshold for selecting a circle from votes.


Missing Data in Machine Learning Classification

Labelling a missing value as -1000 (assuming a non-missing value is much smaller in absolute value) will not cause problems with anything based on decision trees/stumps (bagged/boosted). You might also have some success with mixture models, for example a GMM - the learning algorithm will allocate some components to cover the missing data. Same for non-parametric density estimators, nearest neighbours, or SVM with RBF kernels.

Note that in all these cases, even if the model works, there will be a waste of resources ("model space budget") if feature #25 is always present in your testing data. In other words, your model will have learnt how to deal with tricky examples that will never occur when classifying new data! But if feature #25 will sometimes be missing in some of the data you will encounter "in the field", this is the way to go.

Do not expect things to work with algorithms such as neural networks or support vector machines with linear or low order polynomial kernels.

If feature #25 will never be missing "in the field", it might be simpler to just extrapolate feature 25 in your training data. To do so, you can proceed as follows:

Option 1: Nearest neighbor

Normalize all features (except #25).
For each training example in which feature #25 is missing, search for the nearest training sample, according to features (1 .. 24 and 26 .. 36), for which feature #25 is present.
Assign to feature #25 the value taken from the nearest neighbour.

Works best when you have lots of data (but not too much so as to make the search feasible) with kinky distributions.

Option 2: Regression

Use any regression method of your choice to find a relation ship between feature #25 and the other features, on the subset of your dataset in which feature #25 is not missing. For example, a linear regression.

Use the result to predict feature #25 on the remaining training instances.

Works best when your features are "well behaved" for the regression method you plan to use (for example, see the prerequisites for linear regression).


Resources for learning about filters from scratch?

The best book I know that exclusively treats digital filters is Digital Filter Design by Parks and Burrus.
But any good textbook on signal processing has one or more chapters about properties and design of digital filters. And since the most common IIR filter designs are based on analog filters, these books usually also treat the classic analog filter approximations like Butterworth, Chebyshev, Cauer (elliptic), and Bessel.
Two good undergraduate texts are Digital Signal Processing by Proakis and Manolakis and Digital Signal Processing Using Matlab by Ingle and Proakis.
Also take a look at the answers to this very related question.


How to cut off the initial disturbance automatically?

I just come up with a solution when I am writing question...

I can compare the amplitudes of this period and last steady periods (for my case there is a fixed responding frequency). When the relative error is small enough for 5 or 10 times, it can be considered as a good cutoff point.

below a Matlab script I used, and it worked well for my case.

function [t1,y1,t2,y2]=autocut(f1,t1,y1,t2,y2)
    errlim=10e-2;
    cut1=findcut(f1,t1,y1,errlim);
    cut2=findcut(f1,t2,y2,errlim);
    cutp=min([cut1 cut2]);
    t1=t1(cutp:end);
    t2=t2(cutp:end);
    y1=y1(cutp:end);
    y2=y2(cutp:end);
end

function cutnt=findcut(f1,t,y,errlim)
dt=t(20)-t(19);
T1n=ceil(1/f1/dt)*2;%T1n points in one period, number of periods for A=(max-min)/2
tn=length(t);
nperiod=floor(tn/(T1n));
%% Ast is the amplitude of last 5 periods with steady vibration
tstartn=(nperiod-5-1)*T1n+1;tendn=nperiod*T1n;
Ast=(max(y(tstartn:tendn))-min(y(tstartn:tendn)))*0.5;
%% if the relative error is continuously small for 15 periods,
s=0;%count periods that satisfies the condition, 5 statisfaction periods define the steady
for i1=1:nperiod %
    tstartn=(i1-1)*T1n+1;
    tendn=i1*T1n;
    yamp(i1)=(max(y(tstartn:tendn))-min(y(tstartn:tendn)))*0.5;
    err(i1)=abs((yamp(i1)-Ast)/yamp(i1));%relative error
    if err(i1)<errlim
        s=s+1;
        if s==5
            break
        end
    end
end

cutnt=tendn;


which spatio-temporal properties does dynamic textures exhibits?

Stationarity over time just means that if you come back and look at the traffic sometime later it will look pretty much the same as it did when you first saw it. Sure there will be different cars but even the new cars will look similar to the first ones. The cars might be in different positions but their random arrangement will have the same properties as the original arrangement.


Question about SNR in continuous and discrete domain

In your $snr$ expression, if you put a $1/N$ ratio in front of the numerator summation you'll have the standard expression for SNR with regard to sampled data (digital data).


Representing data in combination of QAM and BPSK

Let's us take the simple case $2\times 2$ MIMO. Let $y_1$, and $y_2$ be the received signal at the first and second receive antennas, respectively. Then we have

$$y_1 = h_{11}x_1 + h_{21}x_2 + n_1$$

and

$$y_2 = h_{12}x_1 + h_{22}x_2 + n_2$$

where $h_{ij}$ is the channel coefficient from transmit antenna $i$ to receive antenna $j$, $x_i$ is the transmitted signal from transmit antenna $i$, and $n_i$ is the AWGN at receive antenna $i$.

We can write the above signals in matrix-vector form as

$$\mathbf{y} = \mathbf{H}\mathbf{x} + \mathbf{n}$$

Apparently, there is interference here. So, we need to do some sort of equalization. Let's use zero forcing (ZF) for simplicity. In this case you get

$$\mathbf{H}^{-1}\mathbf{y} = \mathbf{\tilde{x}}=\left[\begin{matrix} \tilde{x}_1\\\tilde{x}_2\end{matrix}\right]$$

where $$\tilde{x}_i = x_i + \tilde{n}_i$$

for $i = 1,\,2$.

Now, suppose that $x_1$ is 4-QAM, and $x_2$ is BPSK, both have the same symbol time. Will that be a problem? and why? They seem separable to me, and can be detected given that the receiver knows in advance that $x_1$ is drawn from a 4-QAM constellation, while $x_2$ from a BPSK constellation. Actually we do this in OFDM. It's called bit loading, where different subcarriers are modulated with different modulation sizes, and schemes when the CSI is know at the transmitter. Of course in OFDM the symbols are separable because of the orthogonality of the subcarriers, but we can somehow separate the signals in spatial multiplexing, although not perfectly.


The second way is how it is done.
In the fast-time (or range-bin dimension) you are right to perform the matched filter.
However, Doppler information is gathered from sampling pulse-to-pulse. The MTI filter essentially subtracts pulses so that if they have similar phase, or the phase is not changing at all, then there is little to no Doppler being generated and thus the target is not moving.
In other words, the three-pulse canceller applies a filter across the Doppler dimension for each range bin, as the second pulse matrix suggests.


This sounds like the classical job of a polyphase filterbank, and a cascade of polyphase synthesizers. Yes, do this in complex baseband.

First, you use your filterbank to divide the 25 MHz in 25 kHz wide channels. Then, you use polyphase synthesizers to combine subsets these channels back to the channel widths you need, and then you use another one to combine these back to 25 MHz of bandwidth.

Another option would be to simply have a prepared table of filter coefficients for your potential subband widths, and just shift these in frequency (which is rather trivial) accordingly to get out your subbands of choice, then combine these back together with a bog-normal adder.

I'm assuming none of this requires that group delay in all these subbands is identical, for example. Otherwise, the design approach would have to amount for different filter delays.


How do I use a Savitzky Golay filter to find local maxima (in between samples) in a discretely sampled 1D signal?

While I'm not familiar with this specific type of filter, based on the plot you've shown, I would guess that the maxima that aren't found by your process are just butting up against the time resolution inherent in the process. Any kind of "smoothing" implies that there is some time-local smearing of the signal of interest, such that if there are two nearby peaks, it's possible that they will be coalesced into one. It's possible that a lower-order filter might present less of this behavior, likely at the expense of the amount of smoothing that you get.


Time for a preamble (underwater acoustic communication)

In communications, no matter whether it's radio, underwater acoustics or say powerline comms, there's multiple influences on the permissible duration of a symbol.

If you consider your preamble to be kind of a "super symbol" that you use for detection and equalization / phase correction, then you'd typically want to make sure that the channel coherency time is larger than that symbol's length; otherwise, your channel will change while you're still transmitting that preamble, and thus you will lose sensitivity detecting the preamble, or you'll learn about a channel state that simply is already "gone" the moment your preamble finishes.

However, if your receiver does e.g. OOK/power detection, then you can basically not go "too long" with your preamble – the longer it is, the more gain you'll have.

Now, I'm no expert in underwater communications, but all I know indicates that you might be dealing with a heavy multipath channel, and thus, frequency selectivity. Thus, your preamble detector might not only serve the purpose of detecting the start of a transmission, but also use the preamble to estimate the channel and feed the equalizer. In that case, as said above, a preamble that's longer than your channel coherency makes no sense.

Whether 10 s is longer than your channel is coherent is something I can't say. If you're doing communications in a container of liquid that doesn't change, then your coherency is probably "forever" and there's no upper limit. If you're doing communication across a stream with a lot of large debris rushing through quickly, then 10s might be far too long.


Object identification based on both shape and intensity -- i.e., identify the stars and their diffraction spikes at the pixel level

Here is a first stab.  It is essentially a doubly smoothed contour map.


Kernel vs impulse response: what's the difference?

A kernel is a more general concept, but an impulse response is a special case of a kernel. One usage of the term kernel is to describe an integral transform:
$$y(t)=\int_{-\infty}^{\infty}x(\tau)K(\tau,t)d\tau\tag{1}$$
The function $K(\tau,t)$ is called the kernel of the integral transform.
If you compare $(1)$ to the convolution
$$y(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau\tag{2}$$
then you see that the kernel of the convolution $(2)$ is given by
$$K(\tau,t)=h(t-\tau)\tag{3}$$
Note that the term kernel is also used to describe linear transformations in the discrete domain:
$$y[n]=\sum_{m=-\infty}^{\infty}x[m]K[m,n]\tag{4}$$
Discrete-time convolution is again a special case of $(4)$ with $K[m,n]=h[n-m]$.
In signal processing we regularly use linear transforms such as $(1)$ and $(4)$. The most well-known examples apart from convolution are the Fourier transform, the Laplace transform, the $\mathcal{Z}$-transform, and the Hilbert transform (which is in fact just a convolution).


Sampling and reconstruction of signal in Matlab

You are pretty close. Here is a hint: you need to make sure that your sinc pulses are lining up with your samples. You can check this by breaking it down and plotting individually the sinc pulse train that you are getting. Make sure that these line up as shown in your picture and you should be good to go.

%% Sampling and reconstruction demo
clear,clc,close all;

%% Parameters
F = 30;     % frequency of signal [Hz]
Fs = 2*F;   % sampling rate [Hz]
Ts = 1/Fs;  % sampling period [sec]

%% Generate "continuous time" signal and discrete time signal
tc = 0:1e-4:5/F;        % CT axis
xc = cos(2*pi*F*tc);    % CT signal
td = 0:Ts:5/F;          % DT axis
xd = cos(2*pi*F*td);    % DT signal
N = length(td);         % number of samples

%% Reconstruction by using the formula:
% xr(t) = sum over n=0,...,N-1: x(nT)*sin(pi*(t-nT)/T)/(pi*(t-nT)/T)
% Note that sin(pi*(t-nT)/T)/(pi*(t-nT)/T) = sinc((t-nT)/T)
% sinc(x) = sin(pi*x)/(pi*x) according to MATLAB
xr = zeros(size(tc));
sinc_train = zeros(N,length(tc));
for t = 1:length(tc)
    for n = 0:N-1
        sinc_train(n+1,:) = sin(pi*(tc-n*Ts)/Ts)./(pi*(tc-n*Ts)/Ts);
        xr(t) = xr(t) + xd(n+1)*sin(pi*(tc(t)-n*Ts)/Ts)/(pi*(tc(t)-n*Ts)/Ts);
    end
end

%% Plot the results
figure
hold on
grid on
plot(tc,xc)
stem(td,xd)
plot(tc,xr)
xlabel('Time [sec]')
ylabel('Amplitude')

%% Sinc train visualization
figure
hold on
grid on
plot(tc,xd.'.*sinc_train)
stem(td,xd)
xlabel('Time [sec]')
ylabel('Amplitude')


Why are window functions used when calculating spectral estimates, really?

You are right that no real-world signal will have a "convenient closed form" transform. I also find the quoted sentence misleading, and in my opinion it does not motivate the use of windowing. From your question I believe that you do not actually ask about the pros and cons of different window functions, but if I understand you correctly you ask yourself why people cut signals into pieces anyway if they want to analyze their spectrum. The reason is that many real-world signals are simply too long and/or their length is not known. So you don't want to wait till the signal is 'complete', and you also don't want the complexity of processing such a long signal. So you do block-processing, which reduces computational complexity and gives you quickly a spectral estimate. When you cut the signal into blocks you have many parameters that will influence the result of the spectral estimate: the block length, the type of window, the overlap between blocks, etc. This is a complex topic but you'll find a lot of literature about it. If you have access to Oppenheim and Schafer's Discrete-time Signal Processing I would recommend the chapter on Fourier Analysis of Signals Using the Discrete Fourier Transform. It has a lot of valuable information about block processing and windowing.


Fat32's answer is correct and shows a common pitfall.
The reason that you must do this is because recall that the output of the autocorrelation of a signal of length $N$ is $2N - 1$.
You were performing the FFT with the original sample size of $N = 1000$, effectively destroying necessary information to retrieve the autocorrelation of size $2N - 1 = 1999$.


Why the carrier frequency does not appear after AM modulation?

What you did is multiply two sinusoids with each other:

$$z(n)=\sin(n\theta_1)\cdot\sin(n\theta_2)$$

where $\theta_1=2\pi\cdot 40/1000$ and $\theta_2=2\pi\cdot 4/1000$.
The result $z(n)$ can be written as

$$z(n)=\frac{1}{2}\left[\cos(n(\theta_1-\theta_2))-\cos(n(\theta_1+\theta_2))\right]$$

So you get two sinusoids with the sum and the difference of the two original frequencies, i.e. one with $40-4=36Hz$, and the other one with $40+4=44Hz$. This is exactly what you see in your FFT plot. You use a 1000-point FFT, i.e. the frequency at index $i$, ($i=1,2,\ldots,1000)$, is given by

$$f_i=\frac{f_s}{1000}(i-1)=i-1,\quad i=1,\ldots,1000$$

where $f_s=1000Hz$ is the sampling frequency. You get peaks at $i=37$ and $i=45$, which correspond to frequencies $f_{37}=36Hz$ and $f_{45}=44Hz$, exactly as expected.

I think you made a mistake in the mapping of FFT indices to frequencies.

By the way, with standard amplitude modulation you indeed get the carrier in the spectrum because the modulated signal looks like

$$s(n)=[1+mx(n)]\sin(n\theta_0)$$

where $x(n)$ is your message, $m$ is a scaling factor, $\theta_0$ is the carrier frequency, and the term $1$ which is added to the signal results in the carrier wave in addition to the sidebands. You just multiplied two sinusoids without this constant term, that's why you don't see any carrier. This was already correctly pointed out in chirlu's answer.


The Laplacian filter does actually outputs values of 0 on edges, and positive/negative inside image regions. Edges are detected as zero-crossings of the Laplacian.
This makes sense if you remember that the Laplacian is a second-order derivative, and not a first order.

So, the unsharp filter that you describe computes (identity - Laplacian) convoluted with the image, i.e., the image minus something that is zero on edges. Hence, high frequencies (edges) are not attenuated while lower frequencies are weakened.


How to estimate filter coefficients of an all-pole system?

Half of the information about a system is in the phase response, and half of it in the amplitude response. There's no way to reconstruct a filter if you know only either half, unless you strongly restrict the poles mathematically.


evaluating individual terms of 2D DFT (Goertzel?)

The 2D DFT is separable, so you can conceptually compute it by applying 1D DFTs in both directions. For the Goertzel algorithm you can first apply it for a specific column for each row, and then use that column and apply it again to obtain the individual term. Most likely, you can rewrite the Goertzel algorithm to do this in one go, but the computational savings are quite small (for an $N \times N$ image, you would need $N^2 + N$ multiplications using the separable approach, while a 2D version would give $N^2$ multiplications).


How to make the error of the Gardner timing error detector not depend on signal amplitude?

The maximum absolute error term is roughly 2 times the input power for the Gardner TED. Having a fixed (or predictable) input power, through the use of an AGC, means you can scale the error term to adjust the speed at which corrections will be performed.
With M-QAM modulations (for M > 4), the input power will vary depending on the symbol, and thus the maximum error term will also vary. This means that symbols with less energy will move the timing corrections more slowly than symbols with more energy. Given a proper initial synchronization sequence, this should not be an issue when maintaining the lock. Proper care should be taken to make sure the lowest energy symbol can still make the timing offset move if necessary.


As you may know, orthogonality depends on the inner product of your vector space. In your question you state that:

  While sine and cosine are orthogonal functions...

This means that you have probably heard of the "standard" inner product for function spaces:

$$\langle f,g\rangle=\int\limits_{x_1}^{x_2} f(x)g(x) \ \mathrm{d}x$$

If you solve this integral for $f(x)=\cos (x)$ and $g(x)=\sin(x)$ for a single period, the result will be $0$: they are orthogonal.

Sampling these signals, however, is not related to orthogonality or anything. The "vectors" you obtain when you sample a signal are just values put together that make sense to you: they are not strictly vectors, they are just arrays (in programming slang). The fact that we call them vectors in MATLAB or any other programming language can be confusing.

It's a bit tricky, actually, since one could define a vector space of dimension $N$ if you have $N$ samples for each signal, where those arrays would indeed be actual vectors. But those would define different things.

For simplicity, let's suppose we are in vector space $\mathbb{R}^3$ and you have $3$ samples for each signal, and all of them are real-valued. In the first case, a vector (i.e. three numbers put together) would refer to a position in space. In the second one, they refer to three values a signal reaches at three different times. In this example it is easy to spot the difference. If you had $n$ samples, then the notion of "space" would be less intuitive, but the idea still holds.

In a nutshell, two signals are orthogonal if the inner product between them (namely, the integral I wrote above) is $0$, and the vectors/arrays obtained by sampling them tell us nothing about their being orthogonal.


You need Parseval's theorem. For the discrete-time Fourier transform (DTFT) you have the following relation:

$$\sum_{n=-\infty}^{\infty}|x[n]|^2=\int_{-1/2}^{1/2}|X(f)|^2df\tag{1}$$

where $f$ is normalized by the sampling frequency. For the DFT you have

$$\sum_{n=0}^{N-1}|x[n]|^2=\frac{1}{N}\sum_{k=0}^{N-1}|X[k]|^2$$

So due to Parseval's theorem it is always possible to compute a signal's energy in the time domain as well as in the frequency domain.


Procedural wave generation, frequency shift issue

I suspect that the accumulator value is growing very large, which causes the float's resolution to degrade.  To avoid this you need to occasionally reduce the accumulator to manageable levels.  One way to do that is to subtract by 2*pi whenever you exceed 2*pi.


How to decide whether to use AR or MA for smoothing data?

MA (Moving Average) filter is a FIR filter.
In general it is given in this form -
$$ Y(n) = \sum { w_i * X(n-i) } $$

In inherently it is always stable.

AR (Auto regressive) filter is an IIR filter.
In general, it is given that AF filter is in the form
$$ Y(n) = X(n) * w_0 + \sum { w_i * Y(n-i) } $$

If you want to design a filter who's frequency response window of a very accurate shape it is more manageable or easier to achieve through FIR filters but might be difficult with IIR filters. For example if you need to do equripple filters or Linear phase filters - FIR (MA) filters could be straight in design. However, if the exact shape is not so crucial, using AR one could always produce similarly effective filter (in terms of how steeply gain drops after required cut off), with relatively lesser number of taps.

Here is an excerpt from "Digital Signal Processing" Proakis & Manolakis; this is about the comparison of FIR vs. IIR filter design. Page: 724.

  As a general rule, FIR filters are used in applications where there is
  a need for a linear-phase filter. This requirement occurs in many
  applications, especially in telecommunications, where there is a
  requirement to separate (demultiplex) signals such as data that have
  been frequency-division multiplexed, without distorting these signals
  in the process of demultiplexing. of the several methods described for
  designing FIR filters, the frequency sampling design method and the
  optimum Chebyshev approximation method yeild the best results.

  IIR filters are generally used in applications where some phase
  distortion is tolerable. Of the class of IIR filter, elliptic filters
  are the most efficient to implement in the sense that for a given sent
  of specifications, an elliptic filter has a lower order or fewer
  co-efficients than any other IIR filter type. When compared with FIR
  filters, elliptic filters are also considerably more efficient.


How to determine whether a filter is high/low or band pass from the Z transform?

HINT: There's no need to compute the frequency response by evaluating the transfer functions for $z=e^{j\omega}$. Just compute the poles and zeros and see where they are. Estimate the effect of the poles and zeros on the frequency response by trying to see how they influence the behavior of the transfer function on the unit circle.


Convert your input signal to an analytic signal with a Hilbert transformer, if your input signal is real.  Then run a single complex cross correlation filter.  Assuming no frequency offset, the phase of the correlator output at the correlation peak will be the input signal's phase shift from your reference signal.


On the time-frequency relationships of median filters

The median filter is actually an example of a non-linear filtering operation.

This stands in contrast to linear filtering operations (the 'classical' way that involve the convolution of a filter's co-efficients with a signal). Since convolution in time is the same as multiplication in the frequency domain, you can go back and forth between the domains to figure out co-efficients of your filter in either basis.

However, since medial filtering is a non-linear operation, you cannot do this. There is nothing you can convolve your signal with (ie, no linear operation you can do on your signal) such that you would get the same result as you see by using the non-linear medial filter.


You can use the cross correlation to measure the time delay between a signal and a delayed version of itself.

You can use cross correlation to estimate  the transfer function.

You probably would benefit the most from PHAT

https://www.mathworks.com/matlabcentral/fileexchange/26256-calculate-ambiguity-function-and-utilities?focused=5143816&tab=function


How do you interpolate between points in an image (2D), e.g. using splines?

There are various methods on 2d interpolation (this one, and this one). But most of them considered at least 4 points rather than 2. The simplest 2d interpolation is 3 1d interpolation, in which you interpolate the points between (x1-d1, y1-d2) and (x1+d1, y1-d2) as (x2,y2), then you interpolate the points between (x1-d3, y1+d4) and (x1+d3, y1+d4) as (x3,y3). Finally (x1,y1) is obtained by 1d interpolation of (x2,y2) and (x3,y3).


Which anatomical model to use as "ground truth" for validating an MRI segmentation method?

Your question practically answers itself in the "If your method attempts to estimate the partial volume then the fuzzy model is probably what you want. Otherwise, the discrete model will probably do..." quote.

So: if you're only going to segment each voxel into one of white/gray/CSF then yes, the discrete one seems the appropriate one to use.


You need to cross-correlate those two wavfiles : from wikipedia "cross-correlation is a measure of similarity of two series as a function of the displacement of one relative to the other"

The result obtained will give you a measure of the similarity between those two audio signals depending on the time.

Therefore you will need to find the time for which the cross-correlation is at its maximum (ie : maximum similarity), it will indicates you the time delay, phase between those two signals.

A couple of advices :

1/ If your signals are periodic, you only need to compute cross-correlation over one period (or more if you want to average the noise).

2/ Don't forget to take the absolute value of the cross-correlation while searching for the maximum.


Approximating the frequency response of a median filter

For a start, any non-linear system will not have an easily-identifiable frequency response. So, it's really a nonsensical question. I intend no offense; nonsensical questions are often the  most enlightening!

However one way to try to answer your question is to assume that the LTI filter involved is the mean (rather than the median) of the windowed data.

Then your question:

  Under what conditions or for what class of signals might this approximation be ballpark "close"?

becomes:

  Under what conditions or for what class of signals might the mean be ballpark "close" to the median.

In that case, for a purely stochastic signal, the mean and median are similar when the probability density function (PDF) of the signal is symmetric about the mean.

  For what class of signals might this approximation be very inaccurate?

When the signal's PDF is "very" asymmetric.


Why are the Continuous Wavelet Transforms of the same signal drastically different?

There are differences in axis-scales. However, I suspect that the main difference is that the first picture relates to either a real CWT (or time-frequency transform) or the real-part of a complex CWT (or time-frequency transform, again), while the second seems to be a magnitude.
In the first  image, the apparently wiggling part could "planified" by combining it to some imaginary part. Then changing the $y$-scale from linear to logarithmic could make the plots look closer.
Apart from representation issues (scales, magnitudes), with admissible continuous wavelets, there exists a wide variety thereof.  Factors that affect scalograms in practice are also the global shape of wavelets (symmetry and anti-symmetry for instance), their regularity, and the ways it is discretized (number and sampling of voices).


You will have to implement the transfer function in pieces, rather than as a single LTI block.  See http://www.mathworks.com/support/solutions/en/data/1-VPJDU


Generate noise with arbitrary covariance

If you want a Gaussian vector $X$ with $E(X)=0$ and $E(X\,X^T)=C$, then find a decomposition $C=AA^T$, for instance a Cholesky decomposition, compute a vector $U$ with standard normal uncorrelated random numbers as components and set $X=AU$.


Why Would Pre Filtering Measurement Data Affect the Least Squares Estimate?

I'm not sure what's you model is.

Let's say it is something like:

$$ y = H x + n $$

Now, using the Least Squares model is optimal (In the MSE sense) when $ n $ is AWGN (It is the linear optimal estimator if the noise is white).

So unless the noise in your model is colored, no gain by filtering the data before applying the Least Squares method.

Now, what if the noise isn't IID White Noise?

If it is white just with different statistics, then use the Weighted Least Squares.
If it is colored you should apply Whitening filter on the data.

So, if the Whitening filter is given by $ W $ then you have:

$$ z = W y = W H x + W n = A x + v $$

Now you have classic LS model, where your model matrix is $ A = W H $ and your white noise is $ v = W n $.
Now you can just apply the regular Least Squares.


Frequency slows faster when it's decreasing

I think the issue is that you aren't ramping down the frequency from A to B, but the phase.  The phase needs to be integrated from the last previous value.
Plot the phase argument and you'll see what you did wrong.  The phase for source activity 1 should end up parallel but above the phase for source activity 2.
You say that the end frequency is 4 in both cases, but what you really have done is constrict both waveforms to have the same number of cycles.  Since source activity 1 starts faster, it must end slower.


Normalizing a vector with Mel Cepstrum Coefficients, Delta Coefficients and Log Energies for both

You could collect a large collection of voice signals, extract the MFCCs, and compute the mean and standard deviation of each coefficient. Keep these numbers and use them once and for all.

Then, you can normalize the sequence of vectors from an utterance by subtracting the mean and dividing by the standard deviation you have pre-computed (do not recompute the mean and standard deviation on each utterance - since each file utterance will have slightly different statistics).

This might be useful if you use the euclidean distance to compare MFCCs - for example when using DTW. However, more advanced statistical modelling techniques do not need this normalization step - a gaussian model, for example, will actually learn the mean and covariance matrix of the coefficients without having to pre-normalize the data (Euclidian distance on mean/std normalized data is just like Mahalanobis distance with a diagonal covariance matrix).

For sake of accuracy, I should point out that this normalization scheme has a flaw... The last MFCC coefficients are known to be very noisy and less informative - this follows from the energy compaction property of the DCT (the less relevant stuff is in the highest coefficient), or from the view of the DCT as an approximation of a PCA on the log mel spectra. Thus, the mean/standard-deviation scheme I have presented might "blow up" the importance of the last coefficients. Machine learning/feature selection algorithms might be robust to this and would discard them anyway, but the euclidian distance is not. A compromise I have used in past research is to normalize by a power < 1 of the standard deviation ; or to scale down, after normalization, the nth MFCC coefficient by $\alpha^n$ where $\alpha = 0.95$. Other normalization schemes implement "liftering" - by multiplying the sequence of MFCC point-wise by a window function.


Need help deriving the energy of a signal

The factor just comes from the substitution $y=\sqrt{2}\omega$ in the integral defining the error function:

$$\begin{align}E_x=\frac{1}{2\pi}\int_{-\infty}^{\infty}|X(\omega)|^2d\omega&=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-2\omega^2}d\omega\\&=\frac{1}{2\pi}\frac{1}{\sqrt{2}}\int_{-\infty}^{\infty}e^{-(\sqrt{2}\omega)^2}d(\sqrt{2}\omega)\\&=\frac{1}{2\pi}\frac{1}{\sqrt{2}}\int_{-\infty}^{\infty}e^{-y^2}dy\\&=\frac{1}{2\pi}\sqrt{\frac{\pi}{2}}\lim_{x\to\infty}\textrm{erf}(x)\\&=\frac{1}{2\sqrt{2\pi}}\end{align}$$


What is the relationship between the sigma in the Laplacian of Gaussian and the two sigmas in the Difference of Gaussians?

I understand that a Laplacian-of-Gaussian filter can be approximated by a Difference-of-Gaussians filter, and that the ratio of the two sigmas for the latter should be 1:1.6 for the best approximation

In theory, the smaller the ratio between two sigmas, the better the approximation. In practice, you'll get numerical errors at some point, but as long as you're using floating point numbers, smaller values than 1.6 will give you a better approximation.

To illustrate, I've plotted a cross-section of the LoG and DoG for a few values of k in Mathematica:

As you can see, k=1.6 is not an ideal approximation. For example, k=1.1 would give a much closer approximation.

But you usually want to calculate LoG approximations for a range of sigmas. (Otherwise, why bother with the DoG approximation at all? Calculating a single LoG filtered image isn't more expensive than calculating a single DoG filtered image.) So the value of k is usually chosen so that you can calculate a series of gaussian filtered images with sigmas s, sk, sk^2, s*k^3..., and then calculate the differences between adjacent gaussians. So if you choose a smaller k, you'd have to calculate more "layers" of gaussians for the same sigma-range. k=1.6 is a trade-off between wanting a close approximation and not wanting to calculate too many different gaussians.

  However, I'm not sure how the two sigmas in the Difference of Gaussians relates to the sigma for the Laplacian of Gaussian. Is the smaller sigma in the former equal to the sigma of the latter?

From the formulas on the wiki page @Libor linked to, you can see that $t=\sigma ^2$, so the approximate a LoG for some sigma, you need two gaussians with sigmas $\sqrt{\sigma ^2+\text{$\Delta $t}}$ and $\sqrt{\sigma ^2-\text{$\Delta $t}}$ (at least in the limit $\text{$\Delta $t}\to 0$). Or, in terms of k:

$\sigma _{\text{Laplace}}=\sigma \sqrt{\frac{1+k^2}{2}}$


Is it possible to infer the characteristic curve of a recording device only from (many) recordings?

Yes, that's definitely possible. The key would be to find recordings where it's reasonably easy to find an acoustic reference. Great candidates for this would be classical or big band recordings. Let's say you have wax cylinder of Beethoven's 5th, you can try to create a reference by analyzing a dozen or so recent recordings of the same piece. The standard deviation would give you a "margin of error" for the reference and than you can simply divide the spectrum of the wax cylinder by the reference spectrum. Repeat for as many recordings as you can find good references for.

  should have an everage spectrum roughly resembling white noise

No. It's certainly not white: closer to pink but with more high end roll off. Hence you need a known good reference to determine exactly what it is.


Is it correct to subtract two signals acquired over different times and trials to remove the common signal in them?

Will averaging S1 over these 60 acquisitions and S2 over 4 (Figure 1) and then subtracting the averaged signals (S2_avg - S1_avg) give me a better representation of response to Stimulus B or should both the signals be averaged across equal acquisitions (Figure 2).

Ideally, both signals should be averaged over the same amount of repetitions and the ideal here would be the biggest of the two (60). This is for each individual point estimate to "suffer" from approximately the same (and as low as possible) variance.

Potentially, you might gain some more insight if you were to use bootstraping, where you would be forming point estimates of each time instant of your time series and then examine these smaller subsets (potentially through a boxplot) to try and "guess" what is the most likely time series to be formed, give these two datasets. But the fact that averaging is suggested in the first place anyway, means that an assumption is already made for the time instance differences to be coming from a gaussian distribution.


Calculating the magnitude spectrum and phase spectrum

HINT:

Magnitude $M(\omega)\ge 0$ and phase $\phi(\omega)$ are defined by

$$X(\omega)=M(\omega)e^{j\phi(\omega)}\tag{1}$$

Note that $(1)$ is generally complex-valued. In your example, the Fourier transform $X(\omega)$ is clearly real-valued. This restricts the possible values of the phase $\phi(\omega)$. What are those two values of $\phi(\omega)$ for which $X(\omega)$ in $(1)$ is real-valued? How do they relate to the sign of $X(\omega)$? I'm sure you can take it from here.


How does an equalizer prevent a sample and hold effect?

Equalizers do not prevent the sample and hold effect, though if you are talking about compensation filters, see the edit below.

There are two types of reconstruction filters (also known as low-pass filters): analog and digital.  An analog low-pass filter can smooth out the output of a DAC.  A digital low-pass filter cannot, since it comes before the DAC.  What it can do is reduce the effect of the sample and hold by increasing the sample rate (assuming it is an interpolation filter) and then getting rid of the aliases via the low-pass filtering.  The effect of the sample and hold is reduced due to the increased sample rate.

EDIT: As Robert points out, the frequency droop that happens as a result of the zero-order hold can be partially compensated for with digital filters whose frequency domain gain is the inverse of the zero-order hold effect.  While this will not eliminate the sample and hold effect, it will mitigate it.


Adaptive filter for noise cancellation when measuring some input

I'm a little skeptical of your derivation and notation as it abstains from telling whether a Mean-Square or a Least-Squares metric of error is minimized, but rather just works on an instantaneous error at sample $n$. But I believe you eventually would devise an LMS or RLS filter, hence either of the two would happen.

So, based on your notation the adaptive system error is $\hat{x}(n)$, which is an estimate of the true signal $x(n)$,  and you have two correlated noises $x_N$ and $H(x_N)$ which are uncorrelated with the true signal. The standard noise cancelling setup.

Had you used the steepest descend based quadratic MSE surface minimization procedure, your instantaneous gradient on the error surface would be something like:

$$ \nabla J(w) = - E\{ e(n) \bar{x}[n] \} $$

where the expectation is due to stochastic MSE framework. In your notation this would become:

$$ \nabla J(w) = - E\{ \hat{x}(n) \bar{x}_N[n] \} $$

Now, if the optimum condition is reached, then the error is minimized, and we can assume that $\hat{x}(n)$ converged to the true signal $x(n)$ and you may replace it. Remember this is an approximate conditions actually...

$$ \nabla J(w) = - E\{ x(n) \bar{x}_N[n] \} $$

and since your noise source is uncrorrelated with the true signal, this expectation is zero, yielding the zero gradient at the minimum point.

In your case, you don't specifically use expectations, so might be harder to show why the gradient would go to zero. Furthermore, as always, be aware that gradient would only theoretically go to zero. Any practical implementation will deviate from this depending on how unrealistic your assumptions were.


The Fourier Transform is (mostly) symmetric in both directions so the same sampling rules apply.

The minimum frequency resolution is determined by the length in the time domain, the same way that the minimum time resolution is determined by the "length" in the frequency domain.

That creates a bit of theoretical problem: in order for the required frequency resolution to be finite, the length in the time domain has to be finite as well. However, a signal that has finite length in one domain has infinite length in the other. So you can't actually sample a signal in both domains at the same time.

In practice, you need to work around this by replacing "must be zero outside the supported interval" by "must be small enough outside the supported interval". So you accept some residual aliasing in either domain, provided you can make it small enough for your specific application


If the spectrum of your I/Q samples is centered at zero then you'll have to perform either AM or FM demodulation before routing any real-valued audio samples to a sound card. For AM demodulation you'll need to implement a complex-input "envelope detector" which produces a real-valued audio signal riding on a DC bias. (In a few days check the web page:

www.dsprelated.com

for a blog of mine describing different envelope detectors.) That amplitude-biased audio signal then needs to be passed through a "DC blocking" filter (a highpass filter whose cutoff frequency is very close to zero Hz) and the output of that filter goes to your sound card.  From Chapter 13 of my DSP book, the most common DC blocking filters look like the following:

where alpha is a number close to, but less than, one.

Now if your I/Q signal (whose spectrum is centered at zero Hz) is an FM modulated signal there are two common ways to perform FM demodulation to produce a real-valued audio signal sequence.  One way is unpleasant because it requires the computation of arctangents and the other way is unpleasant because it requires you to perform divisions. Oleg, search the Internet for "digital FM demodulation" and you will be able to learn about those two different FM demodulation methods.


To answer your first question, yes, you will need to convert the signal into Z-plane. Bi-linear transformation is one way to achieve the desired result. You can even try Impulse Invariant Transform to convert the analog filter to a digital one.


Assuming you are sampling above the Nyquist rate, then an interpolation kernel or equivalent fractional FIR filter on one or both channels (both to match up any non-zero processing delays) might be an appropriate solution.


How does complex data of an FFT and IQ modulation differ?

If the signal frequency is constant and exactly integer periodic in the FFT width, then the transform basis vector for each FFT result bin is the same as a sinusoidal (complex exponential) IQ down-modulator mixer input that starts with a phase of zero at sample 0.  Thus the atan2() results for a signal within that FFT result bin should be identical.

If you only need one result bin of a DFT, then a complex mixer (complex Goertzel) requires less computation than a full FFT result (assuming equivalent numerical precision, etc.).  If you need O(logN) DFT result bins or more, then the FFT is more efficient.

If the frequency is not known to be exactly integer periodic in the DFT/FFT width, then an FFT might be useful as part of an initial frequency estimator algorithm.


Is the Kalman Filter a Best Linear Unbiased Estimator (BLUE) for Heteroscedastic Noise?

Kalman filter is the best linear estimator regardless of stationarity or Gaussianity. Also in the Gaussian case it does not require stationarity (unlike Wiener filter). In the linear Gaussian case Kalman filter is also a MMSE estimator or the conditional mean.


How to approximate a moving RMS without iterating over N samples for each required output?

There is a way to do it, it is really similar to a moving average

$N$ : number of samples per period

$Z$ : Accumulator

$x[n]$ : current sample

$$Z[n] = Z[n-1]+ x[n]^2 - x[n-N]^2 $$

$$x_{RMS}[n] = \sqrt{\frac{1}{N}Z[n]}$$

Basically, you need a delay line to store the previous $x[n]$ samples or better yet the previous $x[n]^2$ samples and you need an accumulator. Both the accumulator and delay line should be set to 0 prior to starting the moving-RMS.

Basically you need need 1 multiplication, 2 additions, 1 division (by a constant, so it is equivalent to multiply by another constant) and 1 square root operation per sample.

Edit : The number of samples per period should be an integer.  There are ways to do it with a non-integer number of samples per period, but it takes more work.

Edit 2 : This scheme will work only if the accumulator is an integer and if the accumulator never wraps around as Robert Bristow-Johnson pointed out. If the accumulator is floating-point, you might be stuck with some value in the accumulator you cannot get rid of, because floating-point arithmetic is not guaranteed to be associative in all conditions. See R-B-J solution in the comments when using floating-point arithmetic.


I am not sure if this is related, but this cross-product is part of the differential equation, something like $\frac{\textrm{d}M}{\textrm{d}t} = \gamma M \times B$... For a true infinitesimal $\textrm{d}t$, this makes sense. For nummerical considerations, this needn't be true. Think of the following:

neglect relaxation for a moment
$\vec{M}_0$ has a certain magnitude (i.e. length)
when you apply an rf pulse, you expect a rotation of $\vec{M}_0$, and a vector's length is preserved under rotations
when you calculate the rotation from the cross product, with several time-steps $\Delta t$, you probably do something like $\vec{M}_{0,n+1} = \vec{M}_{0,n} + \Delta t\;\gamma\vec{M}_{0,n} \times \vec{B}$ (i.e. you compute the change from the previous magnetization and the rf pulse and add it to the previous magnetization vector)
this operation does not preserve the length of the vector: $\vec{M}_{0,n} \times \vec{B}$ is perpendicular to $\vec{M}_{0,n}$. When you add the vectors, the resulting vector is longer than $\vec{M}_{0,n}$ (Pythagorean theorem)

Hence, if you do not take appropriate countermeasures, your magnetization vector will grow. Maybe this is what you see in your simulation: At the end you have an insanely huge signal that does not belong there. You can easily check this via a simple test: Neglect relaxation and apply a huge flip angle, lets say $90000^\circ$ and monitor your magnetization vector length. Certainly it should remain constant. If it doesn't, there might be a bug in your code.

Also try to reformulate your calculations with rotations matrices and check if the results get better.


Circular Convolution as Cyclic Shift Operator

You may solve it by 3 steps:

Show yourself that a Cyclic Convolution with a vector $ \boldsymbol{e}_{i}^{N} $ is a Cyclic Shift Operator $ {T}_{i - 1} \left( \cdot \right) $. Where $ \boldsymbol{e}_{i}^{N} $ is defined as a vector of length $ N $ which all its elements is zero but the $ i $- th element which is 1.
Decompose $ \phi $ into $ 3 $ vectors of type $ \boldsymbol{e}_{i}^{2 N} $.
Use linearity of the Cyclic Convolution to prove your assertion.


The length of the arrays you get back depends on the mode you use to do the decomposition and the wavelet length, which can vary like you see between db1 and db2:

len(cA) == len(cD) == floor((len(data) + wavelet.dec_len - 1) / 2)

If you switch to periodization mode mode='per', you get a different calculation that may be more consistent for you:

len(cA) == len(cD) == ceil(len(data) / 2)

Regardless, decomposition is generally dyadic, so your dyad width is len(signal)/len(cAn).

Once you go muiltilevel (use level='2', etc.), you'll have small variations in length as well. To account for this, use IDWT to reconstruct.


Lower frequencies is not visible in FFT plot when adding into the original signal

What you're assuming about how FFT works is actually correct for Fourier transform, but – despite its name – FFT is not such a continuous all-$\mathbb{R}$ transformation. It is actually just a DFT, discrete Fourier transform, i.e. both the frequency and time spaces are quantised. In your example, you selected the quantisation in frequency space to be $100\:\mathrm{Hz}$, so there exists no $50\: \mathrm{Hz}$ bin.


Why Is Non Linear Least Squares Method from MATLAB and Alglib Gives Different Results on the Same Data?

When you solve Non Linear Least Squares problem of a non convex cost function the end solution (Which is guaranteed to be a Local Minimum) will depend on:

Method of Minimization.
Method Parameters.
Starting Point.

In the case above you set the starting point to be the same for both.
Yet in AlgLib you use the method of Levenberg Marquardt (Classic for Non Linear Least Squares).
MATLAB used to use Levenberg Marquardt as its default in the past.
Yet in recent versions it uses more modern method called Trust Region.
The trust region based methods limit their step size to be more conservative.
Basically like the LM they approximate the problem using Quadratic Function yet in each step the find the minimum with constrain about the domain (The Trusted Region) hence usually ends in the closest local minima while LM might skip it.
If you share the values of a and b in a single simulation we'll be able to look into the Cost Function (As it depends on 2 parameters).
If you are after more modern methods to solve the problem in C (Which I recommend), you should look for Trust Region solvers:

Gnu Scientific Library (GSL) - Non Linear Least Squares Fitting.
Ceres Solver - Solving Non Linear Least Squares.
OptimLib - Conjugate Gradient (It's not optimized for Non Linear Least Squares, yet can be used).

If you look in GitHub for Non Linear Least Squares and Trust Region you'll find more.
Few tips:

The problem isn't convex hence a good strategy is perturbate the starting point. This is a poor man Global Optimization strategy.
If you're sensitive to the result (Must have the global minimizer) you should start with a rough grid search to get better starting points.
At top level if needed use full Global Optimization framework like you have in MATLAB.


Recommendation for good ressources on digital processing and D-QPSK

Wow, I'm honored by Matt L. doing what I'm often doing: Referring people to GNU Radio.

The project actually has a list of recommended literature, but I don't know how well that'd fit you. It's probably still worth looking into.

Then, regarding QPSK: Well, it's one of the basic constellations, and you'd probably be best off reading a textbook intro to digital communications – I like Proakis' Digital Communication, personally, and the third edition (which sadly reduced the formula count) is available cheaply used from libraries.

Regarding DAB+: It uses a differential QPSK, but not directly on the air, but to drive the subcarriers of an OFDM signal, if I'm not mistaken; so you should understand both, as well as (at least superficially) the channel coding used. This again calls for some textbook knowledge!

There's already a Free and Open Source DAB+ transmitter implementation – opendigitalradio is the home of that.

There used to be a GNU Radio DAB (+) transceiver, maintained by Andreas Müller (not related to me, AFAIK), but I think it's for a long-gone/unsupported version of GNU Radio, so compiling that on a modern system would be a challenge already. However, looking at the source code and flow graphs might be pretty enlightening already! Mr. Michael Höhn implemented he receiver as his Master's thesis, and he came up with this structure (p. 67):

In my head, this diagram pretty neatly translates to a GNU Radio flow graph – that is, a graph built from GNU Radio signal processing blocks. GR is a system to exchange data between these signal processing modules and to implement such yourself. As such, it lends itself very well to implementations of receivers that understand their operation as flow of signal and data through multiple processing stages. In fact, there's a very popular graphical flow graph design tool (the GNU Radio companion) which represents this abstract idea of a signal flow graph as an actual image of connected rectangular blocks. I think, just to get an idea of what SDR implementations might look like, the GNU Radio Guided Tutorials might be a good start.

Furthermore, if you're interested in audio broadcasting, DRM (digital radio mondiale) might be of interest to you, too – and with gr-drm there's a working GNU Radio transmitter implementation for that.

I really do not want to stand between your enthusiasm and diving heads-on into software defined radio and digital signal processing – but you need to be aware that of all engineering disciplines, these are among the more mathematical, algebraically influenced. So your initial "I want to understand the whole machinery of DSP" might be a little too ambitious – I don't think there's anyone on here claiming that they'd do that – it's a huge field of expertise, and everyone can just know the basics and be an expert on their own specialization. There's audio processing, image processing, software defined radio, biosignal, radar, sonar, … experts on here, and some that I look up to as being pretty universally educated in DSP, but none of these is expert in everything. That is a good message for you – you don't have to be an expert in everything to be good in what you do!

So, personally, I might recommend starting slowly, using a Book on Signal Theory, which nowadays always introduces to the basics of DSP, too. Based on that, go for digital communications. You don't have to fully understand every bit of e.g. the information theory used in channel coders and source coders (source coders==compression algorithms), but it's necessary that you understand their role, so having a background in signal theory is not really optional if you want to understand how to recover your data from a very noisy observation.

Have a lot of fun!


How to find out the transfer function of a FIR filter?


frequency response of a time domain subroutine in matlab

The easiest would be to do an FFT of the impulse response. Something like this
%% calculating the transfer function of a time domain LTI model

% create an FFT grid
fftLength = 1024;
% make a unit impulse
diracDelta = zeros(fftLength,1); diracDelta (1) = 1;
% run your "black box" function, to create the impulse response
h = myFunction (diracDelta); % plus whatever arguments you need
% FFT
transferFunction = fft(h);

You can certainly encapsulate this in a function but you need to pass in your black box function as a function handle and need to figure out how to manage additional arguments and state keeping (for filters with state). Given that the code above is so simple, it doesn't seem worth the bother.
Keep in mind that this ONLY works for LTI system. For anything else, the transfer function is undefined.
The choice of the FFT length is important: to avoid time domain aliasing it needs to be long enough to contain the entire impulse response or at least where it's decayed to "low enough".


Aquila DSP C++ Library - Wave file FFT analog frequency off by factor of 4?

Not familiar with Aquila in particular, and I don't see how you're mapping the output spectrum samples to an X-axis that ends at 800 Hz,  but remember that a real-only FFT mirrors its output data in the second half of the buffer.  You may need to run your loop from 0 to window_size/2 if Aquila's FFT is implemented the usual way.

That could account for a factor of 2x.


Normalized cross correlation via FHT - how can I get correlation score?

Ok, I finally found the reason for my confusion - thanks to this answer.

There were a number of reasons why I couldn't get my spatial NCC and frequency domain peak values to agree numerically. Some obvious, some subtle, all stupid. For the sake of completeness, here they are:

I had missed an offset and scale factor that was introduced when
normalising my image data after importing the pixel values in
preparation for FHT.
I was using a predicted maximum value for the
    peak, from my sub-pixel shift calculation instead of the actual
    peak.
I had a low-pass filter in the pipeline, which - for
    performance - I was applying concurrently with the correlation
    operation. So my peak value was scaled by the filter.

Not being mathematically inclined, I was wrongly attempting to "explain" my problem by assuming that using FHT was somehow inherently different to using FFT. After being pointed in the right direction I found the above issues and finally got the results I was expecting.


Proving that FFT (and by extension any correlation based estimator) is unbiased

A had a long email discussion about the topic of bias in frequency formulas, aka estimators with a really sharp Grad student a while back.  It makes me hesitant to speak authoritatively on statistical matters.
I don't like the term estimator as it is too broad.
Formulas come in two basic flavors, either they are exact or they are not.  Exact means mathematically exact.  So, if you feed it the DFT of a noiseless pure tone (real or complex, different formulas for each), the calculation will result in an answer limited only by the precision of the calculation.
Other formulas are approximations.  Approximations are often times cheaper to calculate and good enough for the job.  In heavier noise, the precision advantage of a more precise formula may be rendered useless, even a disadvantage.  Approximation formulas are inherently biased.  Run a set of tones through them and maybe, like a broken clock, it will appear unbiased for certain tones.
Now run add your run-of-the-mill AWGN to the pure tone signal and test the frequency formulas.  Does the Gaussian shape pass through unskewed and still centered?
The (excellent, Thanks Julien) paper I cited in the comment didn't take quite that approach.
I would be quite interested in another, up to date, comparison with emphasis on CRLB calculation (especially on the real tone case) of all my formulas against the known standards (Start with the cited list, I think he missed Macleods.)
Here is another similar page:
http://www.ericjacobsen.org/fe2/fe2.htm
If you do endeavor to do a comparitive analysis, and there no better way to learn the alternatives, and post it somewhere, please cite it in a comment on this answer.
Frequency on bin (which seems to be what you are asking about in terms of a native bin) cases are sort of a different beast than non-whole cycle per frame tones.
The technique cited in my articles can be extended to any arbitrary number of bins.


You could check out my MATLAB function, which estimates the variance of the noise. It is also mentioned in a Mathworks blog post, in comparison with other methods.

Regarding $\sigma$, you mention in the question, simple answer is yes, longer answer is no. It refers to the variance of the noise normalized with image dimensions as follows:

$$\sigma = \frac{\sigma_0\sqrt{0.5\pi}}{6(W-2)(H-2)}$$

$-2$ is due to the convolution kernel, applied previously. However, this is not a final quantity, giving you the absolute noise level. Nevertheless, for the variance, you could use it as a percentage, just like you mentioned.


why is the spectrogram from matlab different from audacity

I see two possible explanations.

Explanation 1:

Audacity is a non-destructive/non-linear audio editor. It primarily manipulate edit decision lists of audio signals, rarely the signals themselves. So when you nudge a block of audio 1.4s to the right, it does not fill the gap with 0 (silence). Instead, it just stores somewhere "this block plays 1.4s from the beginning of the timeline". For example, I have created here a 30s sine wave tone and nudged it at time offset 27s. Notice that the gap between the origin and the waveform is empty rather than showing a flat line indicating a 0 amplitude.

Now, if you apply the "Mix and render" command, all editing and mixing commands will be actually performed on the signal.

Explanation 2:

Digital silence is a sequence of 0. This causes a problem during the computation of the spectrogram because this involves the computation of log (0) to get to the dB scale. Audacity might use a special color in its palette for digital silence; while matlab might add a very low noise.


These are typical references to RF analog components where the three of these terms are used in the same context. They are often passive devices but can be implemented as active devices (requiring external power) as well. In the passive RF device world they are the following:

Combiner: This RF device combines two or more signals to a single output and is further described as a "Power Combiner". It is typically a passive reciprocal device: it can be used in the other direction as a "Power Splitter".  This can be implemented as a resistive combiner/splitter (for example put three 50 ohm resistors in "delta form" to create a resistive two way power combiner/splitter which will have a 50 ohm match when loaded at all ports, 6 dB loss between each port and no isolation. Often they are implemented with transformers at lower RF frequencies (10MHz to 2GHz), or with 1/4 wave transmission lines at higher frequencies (Wilkinson combiner/splitter) with internal isolation resistor(s) such that they have isolation between the input ports.

Mixer: This RF device multiplies the time domain signal presented at two ports and provides the product out. It typically will have an RF, IF and LO port corresponding to "Radio Frequency", "Intermediate Frequency" and "Local Oscillator Frequency". It is often used in frequency translation and phase detection applications as well as many other uses. This is often implemented as a passive device with Schottky diodes in a ring with balanced transformers as a "double-balanced mixer" but has many other forms. As a passive device it is reciprocal (any port can be an input or output) but there are often different frequency ranges for each port.

Coupler: This RF device has an input port, output port, and coupled port. The signal presented to the input port appears at the output port with minimal loss, while some of that signal is "coupled" off and appears at the coupled port.  The 2-way power divider (which is the 2-way power combiner run backwards) is a 3dB coupler and couples off 1/2 the signal power (although that is exactly what it is, we would typically not call it that). More common couplers are -6 dB, -10 dB, and -20 dB. The theoretical loss of a -10 dB coupler from input port to output port would be -0.4 dB, corresponding to the -10 dB signal that got coupled off. Further, the coupled port is isolated from the output port--- a signal presented there, such as a reflection of the coupled signal or otherwise, would be greatly attenuated (~ 30 dB) at the output port but would travel back to the input port. Therefore you can use a coupler in reverse to inject a signal into a system where the signal is fed into the output, the coupled input is fed into the coupled port and the combined signal will come from the input port with the relative coupling ratio. Similar to the power combiner, this can be implemented with transformers for frequencies from 10MHz to 1 GHz or adjacent transmission line structures at higher frequencies (strip-line and micro-strip forms).  There are also 4 port “bidirectional couplers” which have a forward and reversed coupled port. Basically the 3 port devices simply have the reverse coupled port internally terminated, but many of the bidirectional couplers are made as two 3-port couplers back to back for higher isolation performance.


How to Combine Low and High Frequencies of Two Images in MATLAB?

The easiest way doing so would be Laplacian Pyramid.

Yet, it can be done just by using simple Addition operator.
Just add the High Frequency of one image to the Low Frequency of the other.
Keep in mind few things:

Dimensions must be the same.
Otherwise, interpolate to the same dimensions.
It is better to use HPF which is built from the same LPF used.
To the least, they must have the same Cut Off Frequency.
Use floating point arithmetic for this procedure.


Choosing inverse Z-transform equation, given that $|a|<1$

Which one you choose is quite arbitrary: it depends on which region of convergence (ROC) you choose.

As Example 2 and Example 3 of that wikipedia page shows, if you want the signal to be causal, then choose the first one. If you want the signal to be anti-causal, choose the second.

Another way to think about it is what happens when $|z| = 1$. On the unit circle is where the Fourier transform of the signal is defined. If you want that to be finite, then you need to choose the inverse that includes the unit circle.


Why are two cameras mounted in paralle in stereo application?

When you say cameras are parallel  it only means that the center pixels are parallel.

Every pixel has a different angle, and that refers to the overlapping regions you were talking about. Actually there is at most one overlapping voxel (3d pixel kind of) for every 2 pixels!

There is more about it in epipolar geomtry , it's interesting and pretty simple!


I think you gave the reason in your question. A periodic function can be written as a Fourier series, i.e. as a weighted sum of phase shifted sinusoids, where the frequencies of the individual sinusoids are integer multiples of the fundamental frequency. And the odd multiples obviously give you other pitches than the pitch corresponding to the fundamental frequency. This is what causes the same note to sound different when played on different instruments. It's called timbre.

The harmonics of the note C consist of all possible notes (and also in-between notes). The strongest harmonics with pitches different from C are G (the perfect fifth at $3f$), and E (the major third at $5f$).

EDIT: Your comment clarified your question for me. The phenomenon you're addressing is called octave equivalence. The most common argument for explaining octave equivalence is that the perfect octave is the most basic interval after the perfect unison, due to the simple relation between the harmonics of two notes an octave apart. But the assumption of octave equivalence is by no means universally accepted. As an example, every musician knows from experience that certain chord voicings which sound good in some range sound different (usually not as good) in a lower range. For some interesting thoughts on octave equivalence have a look at this page.


Signal distortion under non-linear channel

I think this will help you: Think about how the "frequency" is given by the zero crossings of the waveform, and then ask if the transfer function you show (and for that matter even if it did have sharp corners and clipping), would that have any effect possibly on the zero-crossings? Can we determine the frequency completely and accurately from the zero crossings for a real waveform?
Also consider this: All noise can be decomposed into AM and PM components, meaning Amplitude modulation and Phase Modulation. (Frequency Modulation and Phase Modulation are similar: If you have FM then you have PM since they are related by an integral-- a change in phase over time is Frequency). So if you have a signal that has AM and FM components, and you hard limit it, you remove all the AM components and you would only be left with FM. Here is the thing to think about then: What if you have a FM signal and hard limit it (which is the sharp corners and clipping you mention)-- it would remove AM if it existed, but how would it effect the FM?
For more mathematical thoroughness, explore how ANY pattern that repeats in time over period T can exist in frequency only at frequency n/T for n as the set of all integers. Would the transfer function from x to y that you show possibly change T for any given waveform that repeats over T?


Correct method for up-sampling cross-power spectrum for sub-pixel motion estimation

Preface

I know the OP asked for Matlab, but if you follow our discussion here in the Mathematica chat it seems any answer will be helpful. Therefore, I'm going to present the approach, visualizing each step, so it can be repeated and debugged in Matlab.

Answer

The approach itself is simple and as far as I can see there are only two things to keep in mind. But first things first. I'm using a downscaled version of the Lena image that can be accessed in Mathematica through ExampleData. I shifted it by {-5.3, 8.7} pixels which can be done through ImageTransformation in combination with TranslationTransform.

{dx, dy} = {-5.3, 8.7};

gImg = ImageResize[
   ColorConvert[ExampleData[{"TestImage", "Lena"}], "Grayscale"], 128];
{nx, ny} = ImageDimensions[gImg];
fImg = ImageTransformation[gImg, TranslationTransform[{dx, dy}/{nx, ny}]];

{gData, fData} = ImageData /@ {gImg, fImg};

Keep in mind that shifting the image by fractions of pixel units will resample it and introduce some blur that degrades the image. Therefore, I would not expect that the perfect solution can be found.

What we have in gData and fData are the pixel values of the images which look like this:

Fourier transform

I don't know how fft2 in Matlab exactly works, but usually, taking the Fourier transform of an image does not result where the frequencies are centered. Let me give you an example, by showing how the spectrum looks if you just use FourierTransform

Image[Fourier[gData] // Abs]

To make the upscaling by padding work, we need to center the transform. You can do this, by multiplying each pixel by $(-1)^{x+y}$ which gives a chessboard-like fashion. This is explained in detail in the very easy to read Digital Imageprocessing by Gonzalez and Woods. A pdf of the 2nd edition can be found here and you should follow the discussion on page 154.

In Mathematica, a function to transform an image into this chessboard-like form can be written down as

centerFrequency[data_] := With[{d = Dimensions[data]},
   Table[(-1)^(i + j)*data[[j, i]], {j, Last[d]}, {i,
      First[d]}]
   ];

If we now look at the same spectrum, with this function applied before using the FFT, we get the expected result:

Image[Fourier[centerFrequency@gData] // Abs]

Normalized Cross Power Spectrum

To create an upscaled version of the normalized cross-power spectrum, we follow exactly what is written in the paper

Calculate the centered fourier transform of fData and gData into f and g
Conjugate g
Calculate f*g/Norm[f*g]
Pad the result by a factor k of the original image dimension (I assume a square image in this example)

In code this looks like this

normCrossPow[ref_, tmpl_, k_] := Module[{
  f = Fourier[centerFrequency@ref],
  g = Conjugate[Fourier[centerFrequency@tmpl]], nx
  },
  nx = First[Dimensions[ref]]/2;
  ArrayPad[f*g/Norm[f*g], k*nx]
]

To get the resulting spectrum, we only need to apply InverseFourier, apply centerFrequencies again and look at the real part of the complex array.

spec = Re[
   centerFrequency@InverseFourier[normCrossPow[fData, gData, 3]]
];
Image[spec] // ImageAdjust

Above, I used a factor of k=3.
The white spot in the upper right corner is our maximum we are interested in. Let us find the position of this maximum

max = FirstPosition[spec, Max[Flatten[spec]]]

This gives {36, 22}. Keep in mind, that the y-axis in image coordinates and matrix-coordinates is reversed and that in a matrix the y-values comes first. Additionally, if images have no shift at all, the maximum would be at {1,1} and we need to divide out our factor k again. Note that we divide by 4 (!!) because with a factor of k=1 we would divide by 1, with k=1 by 2, and so on.

({36, 22} - {1, 1})/4.0

This gives {8.75, 5.25} for {dy, dx} but remember that it is actually dy=-8.75 because of the reversed y-axis.

Now, we can use this result and back-transform our fImg

res = ImageTransformation[fImg, TranslationTransform[{5.25, -8.75}/{nx, ny}]]

and look at the combined image of gImg and res

Fits nicely and the only error we see is the black border that was introduced by initially transforming fImg.


Is Bluetooth Low Energy jamming possible with an SDR like the HackRF on GNURadio?

I am not revealing any big secrets here on jamming and anti-jamming techniques, nor would I condone creating any such interference. What I am about to say is quite simplistic and well known, but knowing more details in how jamming can take place and being more educated on it in general can help good actors in minimizing vulnerabilities in future designs.
Yes this is completely feasible but I would do this much more simply:
Create a single sinusoidal tone at 15 MHz and at 39 MHz.
Upconvert this tone using a carrier at 2441 MHz.
This will produce two symmetric sidebands with a tone at:
2402 MHz
2426 MHz
2456 MHz
2480 MHz
So there is one additional channel rather than just the targeted three channels in the hopping approach, but still significantly better and more effective over trying to jam the entire spectrum from 2402 to 2420 MHz, and conveniently the suppressed carrier that is at 2441 MHz is in between channels.
These tones themselves can be modulated if necessary to spread the energy over the channel bandwidth, since single tone jammers can be easily defeated with tone excision. By spreading the bandwidth, it becomes a simple SNR competition and it can't be defeated given sufficient energy, as the duration of the signal occupation will be continuous in time.

If no sources were readily available, below shows a feasible implementation for comparison to alternative approaches. A custom 2 layer circuit board layout would be sufficient and is \$65 from Express PCB , (the sensitive RF areas should be minimized in pcb track lengths, and for this application there would be no issue implementing this design on such a pcb with attention to proper RF layout techniques in those areas)-  so a total cost of \$127 not including the low frequency noise source or DC supplies and DC control voltages. For single tones the low frequency noise source can be replaced with the same circuit as shown for the 51 MHz LO since the crystal oscillator (XO) can also be programmed for a 12 MHz output frequency. The noise source (modulated and filtered 12 MHz waveform) would be ideal for jamming purposes. The optional "carrier nulling" shown is a small DC offset (coupled in with a resistor) that can be adjusted to null any residual carrier feedthrough leakage at 2441 MHz if desired. Not shown in the spectrum above, but this would also result in suppressed carriers from the first upconverter stage in the output at 2413 MHz and 2467 MHz (also conveniently in between channels). The nulling control for this is in having a small DC offset control in the originating noise source. The nulling is feasible (and effective!) since the IF ports of both mixers have a response to DC.


Is it possible to recover a waveform from spectrograms of magnitude *and* phase?

Real/imaginary or modulus/phase are two representations of a complex number that carry the same level of information. Then, a STFT is a redundant mapping from a space of functions over a 1D variable (time) onto a space of functions over a 2D variable (time and frequency). Under mild conditions on the window, there are an infinity of inverses, due to the redundancy. In the continuous setting, the analysis ($a$) and synthesis ($s$) windows only need to overall a little (more pedantic: $0 < \int a(t).s(t)dt <+\infty$). In the discrete form, this can be a little more involved, yet there are several optimized inverses for "almost" every sound redundant analysis filter-bank).

So, if for each point of the $(t,f)$ plane, you have the magnitude in the phase, you can recover the waveform (apart possibly at the extremeties). A side effect of the huge redundancy, inversion even remains (approximately) possible under mild conditions if you only have the modulus or the phase.


What are the implications of resampling the derivative of a signal with a higher frequency than the original signal?

As Laurent said, no amount of interpolation will restore the information that's not in the original $f_s=1\,\text{Hz}$ signal.

It's been discussed in other posts' comments what a suitable digital differentiator looks like; that's not really easy to answer universally. However, it's certainly a high pass filter. So, I'd say: Try to understand your derivate as a filter operation, if possible. Maybe this gives you a suitable approach on how to implement a "good" (by your application's needs) differentiation.

Also, as Laurent said, when you make assumptions on your original signal, you might find a suitable derivative estimator. The typical approach for DSP (and where Nyquist's theorem mentioned above comes form) is to assume that your signal is a finite sum of complex sinusoids, and hence, every digital signal is equivalently representable by its Discrete Fourier Transform. You might find a derivative according to the differentiation theorem of the Fourier transform, if you want to represent your signal as a sum of orthogonal sinusoids, but you might also represent it according to any other base – and derive a differentiation algorithm that fits your needs best.


What you're looking for basically is "Inpainting".

In Photoshop is called "Content Aware Fill".

You should look for articles for Inpainting.

Trivial option would be using "Replication" of the border and later enforce some smoothness (Either by Gaussian Blur or even Poisson Image Editing).

P.S.
We'd be happy if you could join our Community Proposal at http://area51.stackexchange.com/proposals/86832 for Image Processing and Computer Vision.
We need more people to follow, up vote suggested questions with less than 10 up votes and suggest new questions.
Thank You.


How can I plot a 3D graph of a given Laplace Transform of a function?


the problem is asymmetric railing. It only occurs on the Z-axis because it is affected by gravity, when the vibration exceeds the accelerometer's range the asymmetric range causes a bias. When the vehicle is in motion vibration occurs. Solution would be with improved mechanical mounting or artificially clipping the data. See link below:

https://www.analog.com/en/technical-articles/vibration-rectification-in-mems-accelerometers.html


How to estimate for an exponential in WGN

Maximum $d$ you may get is a total number $N$ of $y_n$'s you have - 1. Larger $d$'s yield curvatures that can't be represented by N points.

Solve a matrix equation $A\vec{t}=\vec{y}$ where $\vec{y}$ is a measurements vector (size $N \times 1$) and $A$ is a normalized version of matrix $X$ (size $N \times N$), where each column $A_j$ is calculated as:

$A_{ij} = \frac {X_{ij}} {\sum \limits_j (X_{ij})^2}$

$X = \left [ \begin{matrix}
x_0^0 & \cdots &x_0^{N-1} \\
\vdots & \ddots & \vdots \\
x_{N-1}^0 & \cdots & x_{N-1}^{N-1}
\end{matrix} \right ]$

d you are interested in is a position of largest coefficient in a vector $\vec{t}$ (size $N \times 1$ starting from zero). Your approximation would be $y = t_d \cdot x^d$.

You can achieve better results if your solver finds solution optimal in $L^1$ sense.

This is a sort of exhaustive search, yes.


One simple way for quantification of contract that I can think of is through use of image histogram. Following is my suggestion

Compute Histogram of the Image
From the counts compute entropy

If you just want to try it out you can use the matlab inbuilt function http://www.mathworks.ch/ch/help/images/ref/entropy.html

You can use the entropy value of the histogram as a measure of contract.


Convolution between two vectors. Length and normalization

The scaling difference between your two output is exactly $32767= 2^{15}-1$ which is exactly the maximum amplitude of a 16-bit signed integer. I'm guessing, the difference is how you import the audio into the program: Matlab's $audioread()$ typically normalize the data to $[-1,1]$, i.e. divide by 32767. It would seem that whatever Python method you are using doesn't do this.

Normalization has nothing to do with the length of the filters and it isn't affected by any zero padding. It's all about the scaling conventions of your inputs and outputs, how do you interface with drivers and/or files.


EDIT  The OP mentions below Deve's answer that this was not the "Generalized Cross Correlation" algorithm referred to.  The real one seems to be this one:

This paper seems to give equations about how to implement it in the time domain.

The equation they give is basically what Deve reports below. The key point they make is that the "input" to the cross-correlation is not the bare signals, but the residuals of an LMS algorithm applied to the bare signals.

This as the effect of pre-whitening the input signals.

The method seems to be from this paper:

Azaria, M.; Hertz, D., "Time delay estimation by generalized cross correlation methods," Acoustics, Speech and Signal Processing, IEEE Transactions on , vol.32, no.2, pp.280,285, Apr 1984

doi: 10.1109/TASSP.1984.1164314

Abstract: The problem of estimating time delay by cross correlation methods is reexamined for the whole class of stationary signals. Expressions are derived for the estimation mean square error (MSE) by the cross correlation method, and are shown to be identical to previously published results for Gaussian signals. The generalized cross correlation method is also analyzed, and the optimal weight function for this method is derived. It is shown to be identical to that derived for Gaussian signals by the maximum likelihood method. For the cross correlation method a simplified MSE expression is derived, which is to be used instead of a previously published result.

keywords: {Autocorrelation;Bandwidth;Correlation;Delay effects;Delay estimation;Error analysis;Estimation error;Maximum likelihood estimation;Mean square error methods;Signal to noise ratio}

From which the key formula appears to be:

where


how can I calculate pitch deviation in a large collection of audio files?

Use a very slow LFO to get very long segments of uniform pitch (this will allow you to use long analysis windows and counter 'border' effects), and use the simplest possible waveform for the oscillator to prevent octave errors in pitch detection. Aubio provides an implementation of the YIN algorithm which is fairly accurate.

I personally doubt there's anything fancy going on here - it is very likely that the relationship between the LFO modulation depth and the pitch deviation in cents is linear.


What you've observed is just the difference between baseband modulation (your original plot) and carrier modulation (the second plot). Note that for the AM case, since the signal that modulates the carrier is real, the resulting carrier-modulated signal is symmetric. This means that you can get away with only transmitting one of the sidebands. This approach is called single sideband (SSB) modulation. With an SSB system, the carrier-modulated signal only has the bandwidth of approximately one of the sidebands of the original baseband signal.


How to Remove Gaussian Noise from an Image without Destroying the Edges?

You might need to consider more advanced techniques. Here are two recent papers on edge-preserving denoising:

Edge-Preserving Image Denoising via Optimal Color Space Projection [in color] This paper preserves edges by decomposing the image into an "optimal" color space and performing wavelet shrinkage. The optimal color space belongs to the luminance/color-difference family (think L*a*b*, or YCrCb).
Edge Structure Preserving Image Denoising From the paper:

  Our method is based on [jump regression analysis], and consists of three major steps, outlined below. First, edge pixels are detected in the whole design-space by an edge detector. Second, in a neighborhood of a given pixel, a piecewise-linear curve is estimated from the detected edge pixels by a simple but efficient algorithm, to approximate the underlying edge segment in that neighborhood. Finally, observed image intensities on the same side of the
  estimated edge segment, as the given pixel, are averaged by the local linear kernel smoothing procedure (cf., [35]), for estimating the true image intensity at the given pixel.

(Jump regression models incorporate discontinuities using step functions. The primary author has a book on this subject.)


Which sections to cover for my graduation thesis on blind source separation?

I've taken this course.
I think using these two papers plus those books you mentioned would be so useful in writing a report, thesis , ... .

Blind signal separation: statistical principles

Independent Component Analysis (A Tutorial)

Also you could download my professor's thesis which is available online in this page


Discrepancy in stability conditions when calculating via RH criterion and Nyquist criteria

The full Nyquist criterion can be summarized to the equation $N=Z-P$, where $N$ is equal to the number of encirclements in the clockwise direction of the minus one point, $Z$ is the number of right half plane closed loop poles and $P$ is the number of open loop unstable poles. So for closed loop stability you want $Z=0$, so $P$ counterclockwise encirclements of the minus one point. If the open loop is stable and you want the closed loop the be stable as well, then this simplifies to that you always keep the minus one point on your left, or in other words rotate around the minus one point in counterclockwise direction.

It can be noted that scaling the open loop transfer function by a constant $K$ then the Nyquist plot also gets scaled by this constant. When determining the closed loop stability this is equivalent to scaling the minus one point by $1/K$.


To measure the energy (which is closely related to, but not the same as, "loudness"), calculate the RMS (Root-Mean-Square).

$ E = \sqrt{\frac{\displaystyle\sum\limits_{n=0}^{N-1} s[n]^2}{N}}$

$N$ is the number of samples and $s[n]$ is the sample at time $n$.

You can do this in any block size of samples.  You can do it on all the samples at once to get an average energy, or you can break your sample set into chunks to see how the energy changes with time.

There are also more sophisticated methods that take into account the body and brain's quirks.  You can find more information about that here.


Similar to the determination of the confidence intervals when estimating the population variance from a smaller sample, with a limited amount of data you will have a large variability in your estimates. There really is no way around reducing this error or variability other than taking more data (to the extent it is stationary).
For an optimized approach to measuring the SINAD, I suggest using correlation of your distorted signal with the reference waveform before any noise is added. I detail this approach further at these post:
How can I find SNR, PEAQ, and ODG values by comparing two audios?
Noise detection


Formulating a function on Matlab for the Shannon interpolation formula

Given your samples $x[n]$, here's a code to perform the interpolation by the sinc-functions:

F = 1000; % the sampling frequency used for drawing the continuous function
Fs = 50; % the sampling frequency
T = 1/Fs; % the sampling time

t = 0:1/F:1; % the time-samples for "continuous" plotting
ts = 0:T:1;  % the time at the sample points

% the original function. In your case, this would be integral expression
x = @(t) sin(2*pi*10*t) + cos(2*pi*11*t);
xn = x(ts);  % the discrete sequence

% perform interpolation from xn to x:
interpolated = 0;
for n=1:length(ts)
    % n-1 due to the 1-based indexing of matlab
    interpolated = interpolated + xn(n) * sinc((t-(n-1)*T)/T);
end

hold off;
plot(t, x(t), 'r', 'linewidth', 3); hold on;  % plot continuous
stem(ts, xn);  hold on; % plot sampled
plot(t, interpolated, 'k-o'); % plot interpolated on top of red curve

As a bonus, here's what happens when the signal is not enough bandlimited: I changed the frequency of the first sine from $sin(2\pi 10 t)$ to $sin(2\pi 30 t)$:

Note: The Shannon-Interpolation requires bandlimited signals, so in particular it does not work exactly with time-limeted signals (since these are never band-limited). So, you will see some deviations especially at the edge of the signal.


type dctmtx

you will get:

function c = dctmtx(n)
[cc,rr] = meshgrid(0:n-1);
c = sqrt(2 / n) * cos(pi * (2*cc + 1) .* rr / (2 * n));
c(1,:) = c(1,:) / sqrt(2);

So it's not that hard to generate an non-square DCT matrix, just with minor revision:

function c = dctmtx1(n,m)
[cc,rr] = meshgrid(0:n-1, 0:m-1);
c = sqrt(2 / n) * cos(pi * (2*cc + 1) .* rr / (2 * n));
c(1,:) = c(1,:) / sqrt(2);


An impulse of "area" $A$ means $A\delta(t)$. The nomenclature arises because
$\displaystyle \int_{-\infty}^\infty A\delta(t) \,\mathrm dt = A,$
that is, the "area under the curve" is $A$.

Since the system is linear and
time-invariant, and the response to input $\delta(t)$ is a pulse of height $1$
and duration from $0$ to $T_s$, the response to $A\delta(t)$ is a pulse of
height $A$ and duration from $0$ to $T_s$.  Its "area" is $AT_s$, that is, $A$
multiplied by the area $T_s$ of the response to the unit impulse.