Multiply FS Coefficients by 0 or 1 to get Low Pass Filter?

This is generally better visualized graphically. When you upsample a signal by inserting zeros between samples you get replicas of the original.

Let's take a simple bandpass signal and plot it. All I'm doing is taking noise and lowpass filtering it.

y = filter(fir1(255, 0.2), 1, randn(1, 1024));
plot(20*log10(abs(fft(y))))

Now let's plot an upsampled version. Usually L is used to denote the integer amount of upsampling in some texts, let's go with L = 2 & 4.

subplot(211)
plot(20*log10(abs(fft(upsample(y, 2)))))
subplot(212)
plot(20*log10(abs(fft(upsample(y, 4)))))

Notice how we now have 2 & 4 replicas of the original spectrum. This is effectively what happens when you upsample a signal.

Now the problem is I have spectral content in places that it shouldn't be. This is where lowpass filtering is used. In the time domain I'm filling in the zeros with interpolated data and in the frequency domain I'm trying to reduce that unwanted spectral content to finally arrive at these plots. Notice how the new filter cutoff frequencies are at fc/L since I've effectively increased the sample rate by fs*L. The original bandpass signal went up to 0.2*fs.

subplot(211)
plot(20*log10(abs(fft(filter(fir1(255, 0.10), 1, upsample(y, 2))))))
subplot(212)
plot(20*log10(abs(fft(filter(fir1(255, 0.05), 1, upsample(y, 4))))))


Circular Convolution of Length L of Sequences of Length Greater Than L

Circular convolution is a very close relative to linear one. Both of them are computed based on the same principles. However circular convolution has a periodic result (based on the fact that it's inputs are periodic), hence its computation can be affected by this.

When implementing a circular convolution there is no loss of efficiency if you actually use a linear convolution underhood (vice versa) as their computational efficiencies are very close, except if you are dealing with architecture optimized code.

On the other hand you may be asking for the penalty paid for computing a longer than necessary linear convolution in order to implement an $L$ point circular convolution of $N$ point periodic sequences when $L < N$. In such a case you can improve the efficiency by the following chopping method (borrowed from @DilipSarwate 's a relevant answer ) which I give here a basic MATLAB implementation from which you can see that the linear convolution is performed at the shorter length $L$, instead of the longer one $N$.

Below gives you a way of computing $L$ point circular convolution of $N$ point sequences where $ L < N$.

clc; clear all; close all;

% S0 - Define the Signals
% Note this is meaningful only when L < signal lengths
L = 3;                 % circular convolution length that we want to compute
x = [0 0.5 1 1 0.5];   % signal length is greater than L
h = x;

xL = zeros(1,L);       % these are for short, chopped signals.
hL = zeros(1,L);
for k=1:L              % compute the CHOPPED signals...
    xL(k) = sum( x(k:L:end) );
    hL(k) = sum( h(k:L:end) );
end
yN = conv(xL,hL);     % compute the linear convolution at the chopped length
yL = yN(1:N) + [yN(L+1:end) 0] ; post-process the result.


filter coefficients to know high pass and low pass filter

I think the filters may be sum of a lowpass and a hightpass filter. They may use to extract some high frequency and low frequency components of a signal.
For example to extract the edges of an image and smooth then simultaneously.


L(x,y) here is a function for luma, a non-linear approximation of luminance. As you can see, it's the weighted sum of the R'G'B' values. The prime symbol (') indicates that they have been gamma-corrected for perceptual uniformity. The coefficients used in this case are standardized by ITU-R Recommendation BT.601. Other standards, such as BT.709 (which may be preferred for its higher accuracy), exist as well.

See: http://en.wikipedia.org/wiki/HSL_and_HSV#Lightness. I can't answer whether you can use luma and still say that it's a component of HSL or HSI. If the author of the book gave you that formula, I suspect you can assume so, at least for that text.


What does the Phase part of a FFT of a received signal represent?

The phases of an FFT result are a phase measurement relative to the time that the FFT input sample window was taken.  Each phase measurement is fairly useless unless you know the absolute time that the sample window was taken, or know the offset of the window from another measurement window (another FFT window with a known overlap/offset/delay/etc., or another window sampled at the same time but offset in space/etc., or some combination of the two), or are comparing or using phases across FFT bins (for reconstruction, convolution, etc.)

Whether the phase is positive or negative depends on how the cosine (or evenness or real) and sine (oddness or "imaginary") components of each DFT result bin relate to each other within the FFT aperture or sample window.  Usually a negative phase indicates the opposite oddness component (hi/low versus low/hi asymmetry) than a positive phase.

If you do an FFT-shift before the FFT, then the phase measurement (the evenness/oddness ratio represented by the atan2 function) is relative to the middle of the vector of samples of the window fed to the FFT, independant of whether the signal is exactly integer periodic in aperture or not.


Get the frequency response curve from FIR filter coefficients + sampling rate

As suggested by a comment, scipy.signal.freqz gives the solution. Don't forget to use this to do the rad/sample to Hz conversion.

Here is a ready-to-use code, posted here for future reference:

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
sr = 44100
w, h = signal.freqz(b=[1, -1], a=1)
x = w * sr * 1.0 / (2 * np.pi)
y = 20 * np.log10(abs(h))
plt.figure(figsize=(10,5))
plt.semilogx(x, y)
plt.ylabel('Amplitude [dB]')
plt.xlabel('Frequency [Hz]')
plt.title('Frequency response')
plt.grid(which='both', linestyle='-', color='grey')
plt.xticks([20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000], ["20", "50", "100", "200", "500", "1K", "2K", "5K", "10K", "20K"])
plt.show()

With logarithmic x-axis:

Another example with a comb filter (COEF = [1] + [0] * 200 + [-1], i.e. y[n] = x[n] - x[n-201]):


Relationship between amplitude and sampling rate

In order to resample without aliasing, the resampling process needs to apply an anti aliasing filter at the new Nyquist Frequency (or thereabouts). So chances are, your resampling process applied a low pass filter at 256 Hz or so and that removed A LOT of the energy.


$$x(t)=\begin{cases}
 t+1\qquad & 0\leq t\leq1 \\
 2-t\qquad &  1<t \leq2 \\
0\qquad & \mathrm{elsewhere}
\end{cases}$$

$$f_1(t)=x(t+2)=\begin{cases}
 t+2+1\qquad & 0\leq t+2\leq1 \\
 2-(t+2)\qquad &  1<t+2 \leq2 \\
0\qquad & \mathrm{elsewhere}
\end{cases}$$

$$f_1(t)=x(t+2)=\begin{cases}
 t+3\qquad &-2\leq t\leq-1 \\
 -t\qquad &  -1<t \leq0 \\
0\qquad &\mathrm{elsewhere}
\end{cases}$$

$$f_2(t)=2x(t+1)=\begin{cases}
 2(t+1+1)\qquad & 0\leq t+1\leq1 \\
 2(2-(t+1))\qquad &  1<t+1 \leq2 \\
0\qquad &\mathrm{elsewhere}
\end{cases}$$

$$f_2(t)=2x(t+1)=\begin{cases}
 2t+4\qquad & -1\leq t\leq0\\
 2-2t\qquad &  0<t \leq1 \\
0\qquad &\mathrm{elsewhere}
\end{cases}$$

Now add them
$$y(t)=x(t+2)+2x(t+1)=\begin{cases}
 t+3\qquad & -2\leq t<-1 \\
3t+7\qquad &t=-1\\
 t+4\qquad &  -1<t \leq0 \\
2-2t\qquad &0<t\leq1\\
0\qquad &\mathrm{elsewhere}
\end{cases}$$

Now check $f_1(-1)=2$,$f_2(-1)=2$,also $y(-1)=f_1(-1)+f_2(-1)=4$


There is any way to assign a a "level/histogram bin" to each pixel without iterating through the entire image?

(emphasis mine)
No, because what you describe is an operation that needs to be applied to each pixel.
Also, a histogram needs to be calculated taking all pixels into account (that's the definition of a histogram!).
Anyway, that's not a problem at all – in SAR processing, this operation is definitely of "a lesser concern" when it comes to computational intensity.

Calculate a list of pixel values sorted by intensity
assign value ranges based on that

Sort the pixels
That's inherently an operation that can't be fully parallelized – in the end, there is one list for the whole image. However, you can parallelize it for the most part: Let your multiple processors (GPU threads, …) sort disjunct parts of the image, then merge the sorted lists in order (and that can be done in cascaded ways, too). Pseudocode:
N := numbers of pixels
I := image considered as linear array of pixels
O := output list of sorted pixel values, same length as I
T := temporary list of list of pixels
P := numbers of processing cores

FOR p IN P PARALLEL:
   T[p][:] := sort(I[p*(N/P) … (p+1)*(N/P)])

merge(T[0…(P-1)])

Quantize
Then, just take the sorted list of values, divide into the desired number of bins of equal size, note down the "starting" value of each bin (i.e. the values at O[0], O[N/n_bins], O[2*N/n_bins]...).
Then, (again, in parallel if desired – there's no data interdependence between pixels) move through your image and replace the pixel values by the bin number.


Deep zeros in the spectrum of the input data?

"Deep zeros" in the input spectrum are caused by a channel which strongly attenuates one or several frequencies, often due to multipath fading. This causes problems for a (linear) equalizer because it tries to compensate for these "deep zeros" by strong amplification of the corresponding frequency bands, which leads to noise enhancement. Non-linear equalizers are better at dealing with such problematic channel conditions.


Z-Transform vs. Fourier Transform convergence

There are several signals for which the $\mathcal{Z}$-transform doesn't exist but the (discrete-time) Fourier transform does. One important example is the complex exponential $x[n]=e^{jn\omega_0}$ extending from $-\infty$ to $\infty$. As a consequence, also sinusoidal signals only have a Fourier transform but no $\mathcal{Z}$-transform, as already mentioned in GKH's answer.

Another important group of sequences with only a Fourier transform but no $\mathcal{Z}$-transform are impulse responses of ideal frequency-selective filters, such as ideal low pass, high pass, or band pass filters. The same holds for the impulse responses of the ideal Hilbert transformer and the ideal differentiator.

In general, two-sided signals with a constant envelope (such as the complex exponential and sinusoids), and two-sided signals that decay only as $1/n$ (such as impulse responses of ideal filters) don't have a $\mathcal{Z}$-transform, but they can have a Fourier transform. The corresponding expressions for the Fourier transform may contain Dirac impulses, or they are discontinuous.


Bandwidth of a finite impulse response filter

Yes, the length of a filter limits the "frequency-selectiveness". If you remember modern physics, that's mathematically very related to the Heisenberg uncertainty principle: something that is temporally short can't be very sharp in the frequency domain, and vice versa.

You can find one example of such a constraint in this answer, though the formula only applies to specific filter constraint.

Generally, filter design is one of the oldest application fields of evolutionary / ML algorithms. You'll be delighted to find that even some of the very "classical", non-ML algorithms for bandpass filter design make extensive use of a loss / reward function while optimizing a filter.

If you're doing modern machine learning for filer design, there's a lot of modern literature, but I'd really start by reading up on the Remez exchange (or just Remez) algorithm, to give you a starting direction. Understanding that relatively simple algorithm and its convergence will give you a solid base!

Especially when later writing a report on it, that would be the most basic benchmark for your ML approach that I could think of. Compare yourself to that!


What is the link between the homography computed over 2 images and the homography computed over the same images upside down?

I think when looking at the formula for the homography matrix:
$$
 H_{ba} = R - \frac{tn^T}{d}
$$
where $R$ is the rotation matrix by which $b$ is rotated in relation to $a$; $t$ is the translation vector from $a$ to $b$; $n$ and $d$ are the normal vector of the plane and the distance to the plane respectively (see Homography-3D plane to plane equation).

The relation between the two matrices lies in the normal vector of the plane. So you need to get the normal vector of the plane (out of the homography matrix), and apply the rotation to it, and then compute the homography matrix using the formula above. For the correct decomposing of the homography matrix, you can look at these code samples and this paper.


Why is particular solution zero for an impulse excitation signal?

The correct form of the statement should be

the particular solution is zero for $t>0$

This is simply the case because the input $\delta(t)$ is zero for $t>0$.
So for $t>0$, the impulse response equals the homogeneous solution. For $t<0$, the impulse response is zero because the system is causal. What remains is the moment $t=0$. The only thing that can happen there is an impulse (or its derivatives but let's not get into that here). So the complete form of the impulse response is
$$h(t)=A_0\delta(t)+y_h(t)\tag{1}$$
where $y_h(t)$ is the homogeneous solution of the differential equation. The constant $A_0$ can be determined by looking at the highest derivative of the input on both sides of the differential equation. For the given system we have
$$h'(t)+\frac{h(t)}{RC}=\frac{\delta(t)}{RC}\tag{2}$$
With $h(t)$ given by $(1)$ we see that on the left-hand side of $(2)$ we have a term $A_0\delta'(t)$ with nothing to match on the right-hand side. Consequently, for this specific differential equation we obtain $A_0=0$.
Note that it's easy to find a system for which $A_0\neq 0$. E.g., for a first-order RC highpass filter we have the following differential equation for the impulse response:
$$h'(t)+\frac{h(t)}{RC}=\delta'(t)\tag{3}$$
Using $(1)$ and matching the terms with $\delta'(t)$ on both sides of $(3)$ we obtain $A_0=1$, resulting in
$$h(t)=\delta(t)-\frac{1}{RC}e^{-t/RC},\qquad t\ge 0\tag{4}$$


The average power of a signal $s(t)$ is given by

$$\overline{s^2(t)}=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{-T/2}^{T/2}|s(t)|^2dt\tag{1}$$

If $s(t)$ is periodic, $(1)$ is equivalent to the average power in one period:

$$\overline{s^2(t)}=\frac{1}{T}\int_{-T/2}^{T/2}|s(t)|^2dt\tag{2}$$

where $T$ now denotes the period of $s(t)$.


Detect the beginning of an increasing signal

It is not entirely clear what sort of signal we are dealing with here, apart from the use of the "audio" tag. If the signal had a wider bandwidth then this would be closer to onset detection. But this is not what we are dealing with here. What we are dealing with here is a slow varying waveform that is considered "on", when it has emerged from some kind of background activity. This viewpoint would make it closer to Anomaly Detection or "Outlier Detection" type of problems.

There are a couple of approaches you can take here. One is to fit a model that includes an "activation" parameter and then try to see when does that parameter transits to activation and take that point as the beginning of the onset. If you set off down the path of model fitting, eventually you are going to have to train your model on a number of these curves so that it learns all the different possible "avenues to activation" there might be. For example, the third trace shows background activity, transiting to an intermediate plateau, transiting to full on activation and even in the "activated" region its slope might show further variation.

So, before you start looking at those techniques, maybe you can try the plain old technique of detecting outliers through the statistics of the signal.

As a human being (?) you seem to have a good idea of when you would like to consider this curve as being "on". Therefore, collect all samples of your signal within the "background" region and use something like a boxplot or fit a distribution to this data. The simplest example of that distribution would be a Gaussian that has a mean and standard deviation. This models your "normality" region. Any value that could have emerged from that distribution is dubious as to whether it belongs to the "background" or "activated" segments. But that is not true for all values because soon enough (as time evolves towards the right), the curve will start pushing towards extremal points of the distribution where the probability of generating such a value becomes smaller and smaller.

Putting a hard threshold there, would give you an estimate of where the "activation" region starts.

Hope this helps.

EDIT:

After more information about the problem was shared, I am more inclined to suggest one of the onset detection techniques which will work directly on the audio data.

In any case, the following (cave) illustrations might help a bit more with the earlier suggestions in this post.

The "Human being" comment comes into play when determining the point of earliest transition from "background" to "key on" (rather than doing it automatically). You use the data in the "background" part of the waveform to estimate the statistics of what it means for a sample to be coming from the "background" part and use that to determine a threshold beyond which the samples are now more likely to belong to the "key on" part.

Alternatively:

Combine many takes at similar settings by aligning them on the "Key on" slope and summarise all of this data with a series of boxplots, each one telling you the sort of value limits you can expect at each time instance. Use that information then to choose the point in time when there is deviation from the background.

(Inset images of boxplot and distribution from this and this wikipedia articles respectively.)


Consider a MATLAB/OCTAVE implementation of 1D-DFT/FFT sum:
$$X[k] = \sum_{n=0}^{N-1} {x[n]e^{-j \frac {2\pi}{N} k n} } $$
where $n=0,...,N-1$ and $k=0,...,N-1$

Those identifiers are used to denote the following:

$N$: is the FFT length as in $N$-point FFT.

$x[n]$: is the discrete-time signal whose DFT $X[k]$ is to be computed.

$X[k]$: is the complex-valued DFT of the signal $x[n]$

$n$: is the discrete-time index for the signal $x[n]$ in the range $n=0,...,N-1$.

$k$: is the frequency sample index of $X[k]$, in the range $k=0,...,N-1$.

$j$: is the imaginary unit $\sqrt {-1}$ in signal processing terminology.

Note if the signal $x[n]$'s length is shorter than the FFT length, then zero padding is applied on $x[n]$ to make its length $N$. Otherwise if it's longer than $N$, then the samples after $N-1$ are ignored according to MATLAB/OCTAVE implementation.


Is the Auditory Spectrogram the same as the Magnitude Spectrogram?

Edit: As Jazzmaniac pointed out, the below answer I've given is incorrect. See Jazzmaniac's comment for explanation!

The first source you provide contains this image as an example output of the function audspecgram():

This very much appears to be a standard spectrogram of an acoustic signal. To describe it as an 'auditory spectrogram' might just hint to that fact—that it's derived from an acoustic signal. In other words: you are right in suspecting that what is described as an 'auditory spectrogram' in your source is just a spectrogram that is a result of exactly the procedure you describe to obtain it.


There two definitions for the AM modulated signals
The first one is called as the classical-AM (or conventional AM) and is given by
$$x_{AM}(t) = (A_c + m(t)) \cdot \cos(2\pi f_c t) \tag{1}$$
and the second one is called as the DSB-SC (double side-band suppressed carrier) and is given by:
$$x_{AM}(t) = m(t) \cdot \cos(2\pi f_c t) \tag{2}$$
In your code, you use Eq.2 (DSB-SC) for the first example, and Eq.1 (classical AM) for the second example. Their difference is that the first one includes the carrier spectrum (the middle spike at $f_c$) at the output in addition to the modulating signal spectra (side-band spikes). The latter DSB-SC, however, only includes the modulating signal side-band spectra
and not the carrier spectrum, hence the name suppressed carrier.
DSB-SC has the advantage that the tranmission energy is reduced as the permanently radiating carrier energy is saved yielding a more efficient transmission. Whereas the classical AM, which is less energy efficient to broadcast, has the advantage that its demodulation (by analog hardware) is extremely simple requiring what's known as the envelope detector (a diode, a capacitor and a resistance, in addition ot the antenna receiver and audio amplifier circuitry).


Coefficient weights of Simple Moving Average (SMA) on SMA

The green curve in your plot corresponds to the correct result. The red curve is wrong, because the weights in the vector cnvweights are in the wrong order (left and right halves are interchanged). The correct way to compute those weights is
cnvweights <- convolve(rep(1/n,n), rev(rep(1/n,n)), conj = TRUE, type = "open")
According to the R documentation of the function 'convolve', you need to use conj = TRUE for standard convolution.


Find the stomata in a plant microscopy image

Sorry I don't know OpenCV, and this is more a pre-processing step than a complete answer:

First, you don't want an edge detector.  An edge detector converts transitions (like this dark-to-light):

into ridges (bright lines on dark) like this:

It performs a differentiation, in other words.

But in your images, there is a light shining down from one direction, which shows us the relief of the 3D surface.  We perceive this as lines and edges, because we're used to seeing things in 3D, but they aren't really, which is why edge detectors aren't working, and template matching won't work easily with rotated images (a perfect match at 0 degrees rotation would actually cancel out completely at 180 degrees, because light and dark would line up with each other).

If the height of one of these mazy lines looks like this from the side:

then the brightness function when illuminated from one side will look like this:

This is what you see in your images.  The facing surface becomes brighter and the trailing surface becomes darker.  So you don't want to differentiate.  You need to integrate the image along the direction of illumination, and it will give you the original height map of the surface (approximately).  Then it will be easier to match things, whether through Hough transform or template matching or whatever.

I'm not sure how to automate finding the direction of illumination.  If it's the same for all your images, great.  Otherwise you'd have to find the biggest contrast line and assume the light is perpendicular to it or something.  For my example, I rotated the image manually to what I thought was the right direction, with light coming from the left:

You also need to remove all the low-frequency changes in the image, though, to highlight only the quickly-changing line-like features.  To avoid ringing artifacts, I used 2D Gaussian blur and then subtracted that from the original:

The integration (cumulative sum) can runaway easily, which produces horizontal streaks.  I removed these with another Gaussian high-pass, but only in the horizontal direction this time:

Now the stomata are white ellipses all the way around, instead of white in some places and black in others.

Original:

Integrated:

from pylab import *
import Image
from scipy.ndimage import gaussian_filter, gaussian_filter1d

filename = 'rotated_sample.jpg'
I = Image.open(filename).convert('L')
I = asarray(I)

# Remove DC offset
I = I - average(I)

close('all')
figure()
imshow(I)
gray()
show()
title('Original')

# Remove slowly-varying features
sigma_2d = 2
I = I - gaussian_filter(I, sigma_2d)

figure()
imshow(I)
title('2D filtered with %s' % sigma_2d)

# Integrate
summed = cumsum(I, 1)

# Remove slowly-changing streaks in horizontal direction
sigma_1d = 5
output = summed - gaussian_filter1d(summed, sigma_1d, axis=1)

figure()
imshow(output)
title('1D filtered with %s' % sigma_1d)

The Hough transform can be used to detect ridge ellipses like this, made of "edge pixels", though it's really expensive in computation and memory, and they are not perfect ellipses so it would have to be a bit of a "sloppy" detector.  I've never done it, but there are a lot of Google results for "hough ellipse detection".  I'd say if you detect one ellipse inside the other, within a certain size search space, it should be counted as a stoma.

Also see:

OpenCV: How to detect a ellipse in the binary image
Python and OpenCV. How do I detect all (filled)circles/round objects in an image?
Detection of coins (and fit ellipses) on an image


Is PSOLA or WSOLA better for vocal stretching/retuning?

To be more precise, I'm trying to work out the advantages and
disadvantages of WSOLA vs PSOLA for the following tasks:

Pitch correction:

PSOLA - Advantages you can change the pitch directly based in the rate space between  grains, keep formants intacts...  to change the pitch you not need time scale and  resample after.
PSOLA - Disadvantage do you need a super power pitch detector
WSOLA - Disadvantage WSOLA is not one Pitch Shift algorithm, WSOLA is one Time Scale algorithm, to change the pitch do will need time scale and them resample, the formants will be lost, so you just maybe need one kind of autocorrelation to find nice places to slice, so it is not the same as Pitch information, if are you building an Automatic Pitch Correction one basic autocorrelation will not give to you the correct Pitch Information.
WSOLA - one advantages here is that it is not Pitch dependent to slice, just slice in the best match similarity, remember it is just one time scale algorithm, you will need resample to change the pitch

Time correction:
For time scale I think the WSOLA have best advantage here or some TDHS based algortihm, yeahh you can use PSOLA to time scale and apply pitch shift together, just pitch shift or just time scale , but if are you just time scaling, WSOLA is one smart choice, no need build an Ultra Pitch Detector in WSOLA, just one wave form similariry will solve your problem (AutoCorrelation can solve ...)
Formant correction:

PSOLA - by default formants will be untouchable, but you can warp the formants too, it is amazing in PSOLA ...
WSOLA - after Time scale and resample your formants will be lost and to get the original formants you will need get the formants in the spectral envelope from original signal and apply it in the Pitch Shifted signal, you can do it using an low pass filter in the spectral, use a peak picking and interpolate or use a cepstrum to get the envelope ...

I know that PSOLA is infinitely more complex, but I consider this algorithm really great, if you have a decent PitchTrack it will give you great versatility


General tips for PID tuning of super low-friction actuator?

1 - Do you have some general idea of what your transfer function looks like? Or at least a frequency response? Could you perform some kind of step response? You could then try to fit an order-2 model on your step response and then try to tune your controller.

2 - You mention that your system is underdamped. I assume that you have underdamped stable poles (in the left-hand plane) in your process transfer function. Perhaps you could place zeroes near the poles to damp the system.

3 - I never use a PID with a "real derivative" i.e. derivative = x[n] - x[n-1], the derivative is always filtered in order to stabilize the system. Otherwise modeling errors and high-frequency noise will be amplified by the derivative which is something you don't want. Or instead of a PID you could use a PI with a lead controller cascaded.

4 - I recommend that you consult this book http://www.cds.caltech.edu/~murray/books/AM05/pdf/am08-complete_22Feb09.pdf


Is it possible to spread signal using sum of more than one orthogonal code

The disadvantage is that the signal V21 = V1 + V2 = [2, 0, 2, 0] disappears every other chip interval which makes maintaining carrier and phase synchronization more difficult. Also, while the total energy in V21 is equals to the sum of the energies of V1 and V2, V21 uses 4 times the instantaneous power needed by either V1 or V2 and so the transmitter needs to have greater power-handling capacity which increases the cost. Also, one problem that received a great deal of attention forty years ago was the "near-far" problem which had to do with the fact that different spread-spectrum users received with different power levels caused problems when synchronization was not perfect and the strict orthogonality used in the OP's question was not quite valid; practical stuff that modern engineers don't need to deal with in their MATLAB programs.

Advantages? I will let someone else take on the task of enumerating the advantages of the OP's scheme.


Why does this plot correspond to this function?

The axes in your plot seem to be indexes into an array, as if you were using stem3(Z) instead of stem3(X,Y,Z).

Check the example in the documentation.


Is $ H\left\{x_h(t)\right\} = H\left\{x\left(t-\frac{\pi}{2w_c}\right)\right\} $ where $H$ is Hilbert transform

It is quite unclear what you're actually asking, and most of your equations have no meaning. But one thing that you need to understand is that the equation

$$\mathcal{H}\{x(t)\}=x\left(t-t_0\right),\qquad t_0=\frac{\pi}{2\omega_0}=\frac{1}{4f_0}\tag{1}$$

only makes sense for sinusoidal functions $x(t)$ (of frequency $f_0$). For non-sinusoidal signals $x(t)$ the Hilbert transform of $x(t)$ can only be written as an integral:

$$\mathcal{H}\{x(t)\}=\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{x(\tau)}{t-\tau}d\tau\tag{2}$$

where the improper integral is interpreted as the Cauchy principal value. In general, the Hilbert transform is not just a simple time shift.

For a sinusoidal signal $x(t)=\cos(2\pi f_0t+\phi)$ you can indeed show that the Hilbert transform corresponds to a delay of $t_0=1/(4f_0)$, i.e., in the frequency domain the following holds:

$$-j\,\textrm{sgn}(t)\,X(f)=X(f)e^{-j2\pi f t_0}=X(f)e^{-j\frac{\pi f}{2f_0}}\tag{3}$$


Reconstructing Signal From Its Cyclic Autocorrelation

No, you cannot reconstruct the original signal from the cyclic autocorrelation. The fundamental reason is that it results from an averaging operation. Like the autocorrelation and PSD, the cyclic autocorrelation is a many-to-one functional (you already noted that I think). For example, all BPSK signals with identical timing parameters (symbol clock phase and carrier phase), identical pulse functions, and IID symbols have the exact same cyclic autocorrelation functions. Yet there are an infinite number of such signals obtained by all the different sequences of transmitted bits. In general, there is no single unique signal that gives rise to a particular cyclic autocorrelation.

I note that none of the equations you wrote are the cyclic autocorrelation. In general, it is the infinite time average of the second-order lag product of the signal multiplied by a complex sine wave. Aside from the missing scale factor $1/N$, your third equation is an estimate of the cyclic autocorrelation for discrete-time signals.

There are special cases, though, such as periodic signals. Such signals are non-random. The time-varying autocorrelation function $R_s(t, \tau)$ for such signals is identically equal to the second-order lag product itself $s(t+\tau/2)s^*(t-\tau/2)$. Since the time-varying autocorrelation is composed of the sum of sine-wave-weighted cyclic autocorrelations, one might recover $s(t)$ from the set of cyclic autocorrelations, to within a scale factor of $\pm 1$.

The cyclic autocorrelation is sensitive to phase in the sense that the cyclic autocorrelation for $s(t)$ differs from that for $s(t-D)$ in general. This is a basis for synchronization algorithms.

I'm curious why it is obvious to you that there "should be a tremendous amount of redundancy in the cyclic autocorrelation." Redundant what?


Detecting the rotation of digits in an image

You can use Hough Transform to find dominant lines in the image and then based on rho & theta parameters of the Hough transform, align your text.
First you need to remove unnecessary details from your image through closing as below:

binTextImage=TextImage<30;
close=imclose(binTextImage,strel('disk',20));

Then extract edges,
    edges=edge(close)>0;

Then apply hough transform (the normal form):

H=hough(edges);
imshow(H)

Then find the maximal point in the hough domain,

[X Y]=find(max(H(:))==H)

The rho and theta in the position X,Y represents parameters of the dominant line. Use the theta parameter to find out how much the image is rotated.

More details in MATLAB example in following links is given:
https://www.mathworks.com/help/vision/examples/text-rotation-correction.html


Computer exercises and solutions in signal processing

"Digital Signal Processing: A Computer-Based Approach"  by Sanjit Mitra is what you need I guess, especially the exercises at the end of each chapter. There is a booklet on the Internet again by Mitra, named Digital Signal Processing Laboratory Using MATLAB.
The other option could be Practical Signals Theory with MATLAB Applications.


Which codec and format are Audible downloads?

As this post over at the apple stackexchange says, .aax files from Audible are indeed encrypted m4b files, which means that if I would decode an .aax file without converting it, I would get a .m4b file.
This is possible with the following command:
ffmpeg.exe -activation_bytes XXXXXX -i input.aax -c copy output.m4b

The -c copy ensures that the data will not be converted.


What is the theoretical bit error rate for QPSK Matlab simulation when constellation points are transmitted (no symbol duration)?

TL;DR The theoretical BER of an ideal QPSK system is $\displaystyle Q\left(\sqrt{\frac{2E_b}{N_0}}\right)$ where $E_b$ is the energy per bit and $\frac{N_0}{2}$ is the two-sided power spectral density of the additive white Gaussian noise channel.

The receiver in a standard QPSK system is essentially two BPSK receivers using phase-orthogonal carriers. With a standard matched filter, the signal output in each BPSK receiver has value $E_b$ while the noise variance is $\sigma^2 = E_b N_0/2$.  Putting in the complex-value stuff, the signal output in the QPSK system is a complex number taking on values $(\pm E_b \pm jE_b)$ perturbed by complex Gaussian noise $N_i+jN_Q$ whose real and imaginary parts are independent zero-mean Gaussian random variables with variance $E_b N_0/2$. The decisions that each BPSK receiver makes are independent (because the noise in the receivers are independent), and since each receiver needs to decide between $E_b+N$ and $-E_b+N$, each receiver makes errors with probability
$\displaystyle Q\left(\frac{E_b}{\sqrt{E_bN_0/2}}\right) = Q\left(\sqrt{\frac{2E_b}{N_0}}\right)$. Notice the complete lack of information about the symbol duration $T$ or the channel bandwidth $B$.

"But, but, but..." you splutter, "All this is nonsense.  My receiver
has signal output $(\pm A \pm jA)$ and the noise variance is measured to be $(12+j12)$, not this $E_b$ stuff that you keep on yammering about, and what I want to know is how do I get $E_b$ and $N_0$ when I simulate the system!"

Not to worry. The BPSK receiver is a linear receiver (except for the sampler and decision device), and so the output that you are measuring is really some constant $K$ times $E_b$, and the noise variance $12$ is just
$K^2E_bN_0/2$. So, when you run your simulations, generating gazillions of complex numbers $(\pm A \pm jA)$, perturbing them with complex Gaussian noise of variance $(12+j12)$, your simulated BER should be quite close to $$Q\left(\frac{A}{\sqrt{12}}\right) = Q\left(\frac{KE_b}{\sqrt{K^2E_bN_0/2}}\right)=Q\left(\sqrt{\frac{2E_b}{N_0}}\right).$$
In short, you don't need to know the value of $K$ or to calculate the value of $E_b$ or $N_0$ explicitly; just use your measured values and go from there; it all comes out in the wash.


I see what the upsampling has done but I fail to understand how the amplitude distribution is periodic?

The reason they can say this is because the symbols are assumed to be equally likely. Think about an upsampled BPSK stream, you'll have approximately the same amount of $+1$ pulses as $-1$ pulses so if you take a large enough block of samples it is periodic with period equal to the symbol period. In fact, this periodicity can be exploited to get an estimate of the symbol period (check out the paper Blind symbol rate estimation using autocorrelation and zero crossing detection). For enough symbols, the method does quite well for a blind estimator.

  The author goes on to say "yet within each of the four time offsets the distribution is stationary, leading to the name cyclostationary" : is it stationary within the time offsets?

Yes, as the author says: "within each of the four time offsets the distribution is stationary..." Cyclostationary means that the statistics change for different time offsets but has cycles. Say the peak of your pulse has a value equal to $1.0$. In noise, you sample the peak every symbol (cycle $=$ symbol period here) and you get a bunch of values: $1.1, 0.95, 1.04, 0.98, 0.89$ and so on. You can then collect statistics on these samples like the mean and variance. Cyclostationary means that changing the sample time by a half a symbol period will give you different statistics, but changing the sample time by a multiple of the symbol period (remember cyclo/cyclical) will give the same statistics. This is different from stationary process, where the time shift does not at all affect the statistics.

  "In the frequency domain four perfect images of the original spectrum are created due to the upsampling" how has upsampling caused this to happen?

You can see from the figure you posted, going from the plot labeled "Stationary Symbols" to "Upsampling" we are putting $3$ zeros between every symbol, or every symbol now occupies $4 \times$ as much samples (actual symbol + $3$ zeros = $4$ total samples). This is upsampling by a factor of $4$. When we do this operation, we get a compressed frequency domain axis and $4$ copies of the original spectrum (since we upsampled by $4$). I'm not deriving anything since this is a pretty standard DSP thing, check out
DSP Lecture 14: Continuous-time filtering with digital systems; upsampling and downsampling on Youtube starting at about the 56 minute mark for more in-depth upsampling explanation.


Impulse response of a continuous time system

\begin{align}
y(2.34) &= h \star x\,\big\vert_{t=2.34}\\
&= \int_{-\infty}^\infty h(\tau)x(2.34-\tau)\,\mathrm d\tau
&\scriptstyle{\text{from the definition of the convolution integral}} \\
&= \int_0^\infty e^{-\tau}x(2.34-\tau)\,\mathrm d\tau &\scriptstyle{\text{because } h(\tau)=0\text{ for }\tau < 0}\\
&= \int_0^\infty e^{-\tau}u(1.34-\tau)\,\mathrm d\tau &\scriptstyle{\text{because }x(2.34-\tau) = u((2.34-\tau)-1)}\\
&= \int_0^{1.34} e^{-\tau}\,\mathrm d\tau
&\scriptstyle{\text{because }u(1.34-\tau) = 0 \text{ for }\tau < 1.34}\\
&= 1-e^{-1.34}
\end{align}


What is the maximum value that can result from a 2D DCT?

As suggested in the comments by Marcus Müller, you have to start from the equation of the DCT:
$$
S_{uv} = \frac{1}{4} C_u C_v \sum_{y=0}^7 \sum_{x=0}^7 s_{xy} \cos\frac{(2x+1)u\pi}{16} \cos\frac{(2y+1)v\pi}{16}
$$
The argument of the summation is of the magnitude of $s_{xy}$ times two values which are at most 1 (absolute value). So still 8 bits. You are summing $8 \times 8 = 64$ of those values, so $2^6 \times 2^8 = 2^{14}$. But you also divide this by $\frac{1}{4}$ so the maximum value requires 12 bits. You also have to consider that this value is obtained only when u anv v are 0, and in that case $C_u=C_v=\frac{1}{\sqrt{2}}$, so $C_uC_v=\frac{1}{2}$, that is another bit is unneeded. Some other checks are required when only $u=0$ or $v=0$, but overall you will need 11 bits at most.
This is exactly seen in the JPEG category tables for encoding the DCT coefficients.
For the skeptics, as you can see in the table below, when the input is in [0,255] the maximum values are always below 2048, thus 11 bits are enough. If your input is in [-255,255] you will need 12 bits.

0
1
2
3
4
5
6
7

0
2040
924,25
942,3571
924,25
1020
924,25
942,3571
924,25

1
924,25
837,4883
853,8957
837,4883
924,25
837,4883
853,8957
837,4883

2
942,3571
853,8957
870,6245
853,8957
942,3571
853,8957
870,6245
853,8957

3
924,25
837,4883
853,8957
837,4883
924,25
837,4883
853,8957
837,4883

4
1020
924,25
942,3571
924,25
1020
924,25
942,3571
924,25

5
924,25
837,4883
853,8957
837,4883
924,25
837,4883
853,8957
837,4883

6
942,3571
853,8957
870,6245
853,8957
942,3571
853,8957
870,6245
853,8957

7
924,25
837,4883
853,8957
837,4883
924,25
837,4883
853,8957
837,4883


Make a signal that fits another the best possible with a limitation in the 2nd derivative

Hmmmmmmmmm, interesting question.

Since you want to use the second derivative as your criteria, it would seem that you would want to have the maximum second derivative absolutie value for as short of a duration as possible.  This would suggest piecing together parabolas, matching the first derivatives at the joints.

How to do this algorithmically will take a little more thought.  Somehow you have to identify regions of "violations" and then figure out how big of an interval you have to use to patch it piecewise with parabolas.

Just some thoughts.  If I think of a sure fire way, I'll add it on.

Ced

Followup:

It turns out to be really simple to make a FIR that reacts to a step function with a upward parabola for the first half and a downward parabola for the second half.  It just becomes a matter of selecting the width and the rise rate to set the value of the second derivative.  I haven't done the math yet for the actual value, but here is a program which gives proof of concept.

 h[k] = [ 1  2  3  4  3  2  1 ]  / 16

It appears the second derivative proxy $ y[n-1] - 2 * y[n] + y[n+1] $ has an absolute value of the normalizing factor of the integer step vector, e.g. $1/16$.

import numpy as np

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

        P = 2

#---- Build the Weights

        K = 2 * P - 1

        h = np.zeros( K )

        for k in range( P ):
            h[k] = k + 1

        for k in range( P - 1 ):
            h[P+k] = P - k - 1

        theRecip = 1.0 / sum( h )

        print h, sum( h )

        h *= theRecip

        print h

#---- Set up the Step Function

        x = np.zeros( 20 )

        x[10:20] = 1.0

        print x

#---- Apply Filter

        y = np.zeros( 20 )

        for n in range( 8, 20 ):
            Sum = 0.0
            for k in range( K ):
                Sum += x[n-k] * h[k]
            y[n] = Sum

#---- Show Result

        for n in range( 6, 18 ):
            second = y[n-1] - 2.0 * y[n] + y[n+1]
            print n, y[n], second

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

Yet a better answer:

[Moved to a separate answer]


Inverse Discrete-Time Fourier Transform of $X(Ω)=jΩ$

This system is a discrete-time differentiator, and in this case it's easier if you directly solve the inverse DTFT integral:

$$x[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(\omega)e^{jn\omega}d\omega=\frac{j}{2\pi}\int_{-\pi}^{\pi}\omega\, e^{jn\omega}d\omega\tag{1}$$

Using integration by parts it's quite straightforward to solve $(1)$. The result should be

$$x[n]=\begin{cases}0,\quad n=0\\\displaystyle\frac{(-1)^n}{n},\quad\text{otherwise}\end{cases}$$


Time-spread multiple echoing of audio signals

I'm not that familiar with the technique but reading up it looks like you basically just want an echoic convolution kernel of some sort which you can apply. I don't think the kernel you are building there is correct - if you want, for example, a single echo then the convolution kernel would need to contain two impulses, one at zero offset and one offset into the kernel by the time delay of the echo. For further echoes you would want further impulses, again offset into the kernel by their delay amount. You can have more sophisticated echo kernels but I would suggest this as a starting point

In response to question 2 I would not crop your signal unless it is particularly important to keep the length the same. Keeping the tail will make the deconvolution more accurate and unless your kernal is huge it won't be a very noticeable extension of the signal.


How to choose good features for my machine vision application?

The choice of the algorithm is dependent on the constraints of the application. If you are able to collect (or find) thousands of labeled images covering many possible scenarios, then the popular approach would be to train a deep neural network for semantic segmentation (sidewalk, road, car and etc.). The features are learnt and are not hand-crafted. There are many packages such as Tensorflow or Torch for doing this, as well as trained models. This is one of the technologies powering autonomous driving these days.

If you would need a faster solution and do not have the resources to invest in such data collection, then I would say it is still doable because roads and sidewalks can look very different (depends which sidewalks of course). As the texture of the two classes will be quite dissimilar, I would say that any classical algorithm, capable of capturing texture features could do a decent job. If we are thinking of using OpenCV, then I would immediately try with a HoG (histogram of oriented gradients) + SVM (support vector machines) pair to extract features and classify, respectively. There is already an OpenCV demo (link 1, link 2) for performing object detection in a similar way. Similarly, one could do a sliding HoG extraction + binary classification in order to do a segmentation.

Of course one needs to gather some amount of data for both approaches, but I would say much less for the second one. Rough estimation is difficult, as I haven't seen your sidewalks or roads and cannot predict the variation in the data.


Spreading symbols with sequence before transmission

It is just the control channel that is spread.  It appears that it is spread primarily to make the signal more robust in the face of frequency selective fading, which is probably a major problem at 60 GHz.  It is easy to use spreading to get the robustness because the data rate of the control channel is low.  In other words, the spreading is fulfilling two roles: converting the low control channel data rate to the high OFDM symbol rate, and increasing the control channel's robustness.


How to find the impulse response of filter that outputs the autocorrelation of the input?

HINT:

Note that

$$\int_{-\infty}^{\infty}x(\tau)x(t+\tau)d\tau=\int_{-\infty}^{\infty}x(\tau-t)x(\tau)d\tau$$


Inverse Fourier of Two-Pole Transfer Function

The comments to the question already hint at the solution, but I'll write it out.

$\tilde{H}(\omega)$ can converted from a product into a sum using Partial Fraction Expansion:

$$\tilde{H}(\omega) = \frac{1}{1 + i \omega \tau_1} \frac{1}{1 + i \omega \tau_2} = \frac{A}{1 + i \omega \tau_1} +\frac{B}{1 + i \omega \tau_2}$$

$$ 1 \omega ^0 + 0 \omega ^1 = A(1 + i \omega \tau_2) + B(1 + i \omega \tau_1)$$

Which gives
$$\implies A + B = 1$$
and
\begin{equation}
\begin{split}
A i \omega \tau_2 + B + i \omega \tau_1 &= 0\\
(1 - B) \tau_2 + B \tau_1 &= 0\\
B(\tau_1 - \tau_2) &= -\tau_2\\
B = \frac{\tau_2}{\tau_2 - \tau_1}
\end{split}
\end{equation}

$$A = 1 - B = \frac{- \tau_1}{\tau_2 - \tau_1}$$

Leading to:

$$\tilde{H}(\omega) = \Big(\frac{1}{\tau_2 - \tau_1}\Big) \Big[ \frac{\tau_2}{1 + i \omega \tau_2} -\frac{\tau_1}{1 + i \omega \tau_1}\Big]$$

$\tilde{H}(\omega)$ also contains a term $e^{-i\omega t_d}$, so now for each of the terms in the sum above the inverse Fourier transform is:

$$\mathcal{F}^{-1}\Big\{\Big(\frac{\tau}{1 + i \omega \tau}\Big) e^{-i\omega t_d}\Big\} = \Theta(t) ~ e^{\displaystyle -\frac{t - t_d}{\tau}} $$

and their sum has the desired form in the question.


12kHz sine or square. I can hear the difference - why?

Try synthesizing your “square wave” using only sinusoidal components below half the sample rate.  Otherwise, most of the infinite train of odd harmonics can alias to some other (lower) frequency in the sampled data.


That was the same topic of my Graduation Project, I am not yet sure of the best way to do it as this matter is still a research matter, but what my team and I did was as follow:

1- Provide a database representing the features of the right pronunciation of each word. This database contains for each word around 10 different features for 10 different people reading this word.

2- Get the features of the input word ie. a .wav file representing the word

3- compare the features of the input word with the 10 features representing the right pronunciation of the word

To get the features, you will need to get the MFCC file a .wav file containing the word pronounced
The MFCC is a matrix of numbers, so to compare 2 MFCCs, you will need to use dynamic programming algorithm, in our case we've use Dynamic Time Wrapping algorithm.
You might also need to make a filtering on the input words to filter them from noise.
You might also need to check for the following programs:
Praat  >> to make a filtering if needed, and also can be used to make MFCC file
SPPAS  >> to get the feature of the .wav file (MFCC)
Matlab >> to compare 2 MFCC files

Hope this can help


Determining false positives in feature matching

A simple way to reduce the number of erroneous matches is to tweak matching parameters, including trying different matching strategies. If you are using a simple distance threshold to determine what is a match, you can tighten that. Alternatively, you can use the ratio test (used in David Lowe's demo SIFT implementation) to eliminate ambiguous matches. There you can tightened both the distance threshold and the ratio threshold.

Unfortunately, the world is not perfect, so even if you tweak every parameter, you will still have some false matches. Then you need to use additional information and/or make additional assumptions about your problem. For example, you can assume that the images are related by a rigid transformation (e. g. affine). You can estimate the parameters of this transformation from your matches using an algorithm like RANSAC. Then you can reject the matches that do not fit this transformation as erroneous.

Here is an example of removing erroneous matches by assuming a similarity transformation, using the Computer Vision System Toolbox for MATLAB.


The window length should be equal to your transform length, not necessarily the length of your entire data set.  The two are the same, of course, if you are going to transform the entire data set at once, but if you are planning to do shorter transforms then you should make the window length equal to the length of those transforms.


What is the physical significance of the Fourier coefficients in the DFT of audio signals and how can they be best displayed in a spectrogram?

This is a two year old question, but I think I add a little bit for those who might find it doing searches.

When you are talking about a spectrogram, there are really two steps involved that have to be considered.  First is converting from the signal to DFT bin values, and the second is displaying those bin values in a pixel display.

Ideally, a spectogram of a pure tone of increasing frequency should make a solid, non varying, upward sloping curve.  Inherently DFTs don't work this way.  When you are at a whole integer frequency (in terms of cycles per frame), there is a single bin value, the rest are zero.  When you are halfway between integer frequencies, the bin values are spread out which is commonly known as leakage.  My recommendation is to use a VonHann window on your signal before you do the DFT.  Depending on your scaling, this is the equivalent of subtracting the average of the two neighbor bin values from each bin value.  This has the effect of "tightening up" the leakage, and "spreading out" the sharp peaks.  This makes the display look a lot more consistent independent of the value of the frequency compared to the frame length, which is what you want.

The phase values in the DFT are framing dependent, and therefore don't help much with the display.  Instead you should concentrate on making an intensity display based on the magnitude of the bins.  You are not limited to a linear mapping, or a log mapping.  You are going to want to use a lookup table anyway to speed up processing, so you can pretty much define an monotonically increasing function that works the way you want.  Use can use a monochromatic color scale, or changing color one.  The standard green for low, yellow for mid level, and red for intense is easily comprehended.

The second choice you need to make is the vertical axis scale.  The two logical choices are linear, which the DFT is inherently. or logarithmic corresponding to musical scales.  In either case, you are going to need to map a bin, or a set of bins, to a pixel, or a set of pixels, depending on the magnification at that part of the scale.  Again, in practice, you will build lookup tables for this mapping.

The most proper choices in the above considerations are the ones that give you the visual results you desire.  You are not going to be making calculations based on the display values so there is not problem with "distorting" the values (i.e. non-linear mappings) to give better visual results.

There actually is a physical significance of the values of the DFT coefficients.  Discovering this was a major aha moment in my understanding of the DFT.  A DFT bin value is actually a weighted average of a set of Roots of Unity, where the signal values are the weights and they are "wrapped" around the unit circle.  The bin index is the step size in roots, and thus also then number of times the function is stretched and wrapped around the circle.  For instance, suppose you have a signal with five bumps in it.  When is stretched five times and wrapped around the circle five times, all the bumps align and throw the center of mass in that direction.  On the other hand, when it is strectched by a factor of four, and wrapped around four times, the bumps are evenly spread and cancel each other out.  This is how a DFT "plucks" out the frequencies that are in a signal.

I have written a blog article "DFT Graphical Interpretation: Centroids of Weighted Roots of Unity" to explain this.

This is long enough, hope it helps some one.

Ced


Impulse response with sweep from measurement

referenceFilename = 'reference.wav';
responseFilename = 'response.wav';

[reference, referenceFs] = audioread(referenceFilename);
[nReferenceSamples, nReferenceChannels] = size(reference);

[response, responseFs] = audioread(responseFilename);
[nResponseSamples, nResponseChannels] = size(response);

if referenceFs ~= responseFs
  error('reference audio and response audio files must have same sample rate');
end

nFFT = 2^( ceil( log2(nReferenceSamples+nResponseSamples) ) );

x = zeros(nFFT, 1);
y = zeros(nFFT, 1);

x(1:nReferenceSamples) = reference(:,1);    % use only left channel
y(1:nResponseSamples) = response(:,1);      % use only left channel

h = ifft( fft(y) ./ fft(x) );

h = h(1:nFFT/2);                            % truncate latter half

plot(h);

sound(h, responseFs);


Why Does the DFT Assume the Transformed Signal Is Periodic?

There are already some good answers, but I still feel like adding yet another explanation, because I consider this topic extremely important for the understanding of many aspects of digital signal processing.

First of all it is important to understand that the DFT does not 'assume' periodicity of the signal to be transformed. The DFT is simply applied to a finite signal of length $N$ and the corresponding DFT coefficients are defined by

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N},\quad k=0,1,\ldots,N-1\tag{1}$$

From (1) it is obvious that only samples of $x[n]$ in the interval $[0,N-1]$ are considered, so no periodicity is assumed. On the other hand, the coefficients $X[k]$ can be interpreted as Fourier coefficients of the periodic continuation of the signal $x[n]$. This can be seen from the inverse transform

$$x[n]=\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{2}$$

which computes $x[n]$ correctly in the interval $[0,N-1]$, but it also computes its periodic continuation outside this interval because the right-hand side of (2) is periodic with period $N$. This property is inherent in the definition of the DFT, but it need not bother us because normally we're only interested in the interval $[0,N-1]$.

Considering the DTFT of $x[n]$

$$X(\omega)=\sum_{n=-\infty}^{\infty}x[n]e^{-jn\omega}\tag{3}$$

we can see by comparing (3) with (1), that if $x[n]$ is a finite sequence in the interval $[0,N-1]$, the DFT coefficients $X[k]$ are samples of the DTFT $X(\omega)$:

$$X[k]=X(2\pi k/N)\tag{4}$$

So one use of the DFT (but certainly not the only one) is to compute samples of the DTFT. But this only works if the signal to be analyzed is of finite length. Usually this finite length signal is constructed by windowing a longer signal. And it is this windowing which causes spectral leakage.

As a last remark, note that the DTFT of the periodic continuation $\tilde{x}[n]$ of the finite sequence $x[n]$ can be expressed in terms of the DFT coefficients of $x[n]$:

$$\tilde{x}[n]=\sum_{k=-\infty}^{\infty}x[n-kN]\tag{5}$$
$$\tilde{X}(\omega)=\frac{2\pi}{N}\sum_{k=-\infty}^{\infty}X[k]\delta(\omega-2\pi k/N)\tag{6}$$

EDIT: The fact that $\tilde{x}[n]$ and $\tilde{X}(\omega)$ given above are a DTFT transform pair can be shown as follows. First note that the DTFT of a discrete time impulse comb is a Dirac comb:

$$\sum_{k=-\infty}^{\infty}\delta[n-kN]\Longleftrightarrow\frac{2\pi}{N}\sum_{k=-\infty}^{\infty}\delta(\omega-2\pi k/N)\tag{7}$$

The sequence $\tilde{x}[n]$ can be written as the convolution of $x[n]$ with an impulse comb:

$$\tilde{x}[n]=x[n]\star \sum_{k=-\infty}^{\infty}\delta[n-kN]\tag{8}$$

Since convolution corresponds to multiplication in the DTFT domain, the DTFT $\tilde{X}(\omega)$ of $\tilde{x}[n]$ is given by the multiplication of $X(\omega)$ with a Dirac comb:

$$\begin{align}\tilde{X}(\omega)&=X(\omega)\cdot\frac{2\pi}{N}\sum_{k=-\infty}^{\infty}\delta(\omega-2\pi k/N)\\&=\frac{2\pi}{N}\sum_{k=-\infty}^{\infty}X(2\pi k/N)\delta(\omega-2\pi k/N)\end{align}\tag{9}$$

Combining $(9)$ with $(4)$ establishes the result $(6)$.


How to optimally normalize filter coefficients?

The OP should first implement the filter using floating point to confirm all specifications and structure. Once confirmed there is no "optimal" since the solution is a trade between precision required and accuracy needed in the result but typically the design approach is to use 2 more bits of precision for the coefficients over what is used in the datapath; scale all coefficients accordingly, let the filter grow the signal with an extended precision accumulator (such that no clipping occurs with input set to the maximum possible signal) and then scale the output back to the datapath precision.
See my further details on this at this answer recently posted:
What is the correct gain of an RRC Filter?


Why does Fourier zero-fill interpolation pad the middle frequencies?

Usually a fast Fourier transform (FFT) returns a vector that starts with the zero-frequency bin. This means that the Nyquist frequency (sampling frequency / 2) bin is in the middle of the vector. So you are not pushing apart the low and high frequencies by zero padding at the middle. The zeros go between the highest positive frequency bin and the highest-absolute-frequency negative frequency bin.

Rick omits the details of how to deal with the Nyquist frequency bin (which is there with even-length DFT): If you have a non-zero Nyquist frequency bin then you need to split it into a half-amplitude positive Nyquist frequency bin and a half-amplitude negative Nyquist frequency bin, and insert the zero padding between those. ...I think.


Hint: define a bound for $|x(t)|$, i.e., $|x(t)|\le A$; now try to find a positive number $B$ such that $|y(t)|\le B$ for any $|x(t)|\le A$ (that's simply the definition of BIBO stability). For the system $y(t)=1/x(t)$ it should be easy to show that the above cannot be satisfied (because $|x(t)|$ can get arbitrarily small).


Independent variable transformations (such as exemplified in your question) play an important role in mathematical formulation of signal processing concepts and therefore should be understood well and practiced without difficulties before any real mastery of subsequent concepts. Fortunately they're easy to master! as long as a few tricks has been followed carefully. The most important rule is to break a composite transform into basic transform steps, one at a time.

Method-1:  given $x[n]$ to find $x[3-n]$ proceed as follows:

Define $x_1[n] = x[-n]$ , (reversal of $x[n]$ about origin)
Define $x_2[n] = x_1[n-3] = x[-n]|_{n=n-3} = x[-(n-3)] = x[3-n]$ , (right shift of $x_1[n]$ by 3 samples)

Hence $x[3-n]$ can be obtained from $x[n]$ by first reversing $x[n]$ about origin and then shifting the reversed signal $x_1[n]$ to right by 3 samples.

NOTE-1: While decomposing a composite transform such as $n=3-n$ into its simple steps, you can define a new signal at each step such as $x_1[.]$ and $x_2[.]$, which results in a more cumbersome notation but eventually pays off so that you can clearly see the result of every simple transformation independently.

Method-2: The following set of steps could also be used to get the same result: Given $x[n]$ to find $x[3-n]$ proceed as follows:

Define $x_1[n] = x[n+3]$ , meaning a left sihft of $x[n]$ by 3 samples.
Define $x_2[n] = x_1[-n] \rightarrow x_2[n]= x[3-n]$ , meaning a reversal of $x_1[n]$ by about origin.

Hence $x[3-n]$ can be obtained from $x[n]$ by first shifting $x[n]$ left by 3 samples and then reversing the resulting signal $x_1[n]$ about the origin.

NOTE-2: Be careful to observe the following during step-2: $$x_1[n]|_{n=-n}=x_1[-n] = x[n+3] |_{n=-n} = x[(-n) + 3] = x[-n+3] = x[3-n]$$ and not $$x_1[-n] = x[-(n+3) ],$$ a commonly practiced mistake...


How to analyze multiple amplitude modulated signals?

Well, you can of course instead of using the Discrete Fourier Transform as filter bank just employ bandpass filters to extract your individual tones.
(The DFT is really just a filter bank, if you think about it: it correlates an input sample vector with different complex oscillations. Which is convolution with their conjugate time inverse, which is a filter operation.)
This might be advantageous, since the filter shape of a DFT is inherently sinc-shaped, which I guess isn't quite close to what you want in terms of spectral sidelobes.
So, yeah, simply design a narrow low-pass filter, and shift it in frequency to the appropriate position. (Do that for every tone.).
If you've got many tones that you space on a regular grid, you might want to do a polyphase filter bank, but I'd consider this an "advanced computational optimization that you probably don't need in your measurement process, as you don't care whether computations take a milli- or 100 microsecond".


Fourier transform of a normalized vector

There are different conventions in scaling of the FFT, in MATLAB you need to scale it by $\sqrt{N}$, where $N$ is your number of samples. Saying it in matlabish:

clc, clear all
%% Create some sinusoidal signal with noise
t = linspace(0, 6*pi, 1000)+randn(1,1000);
s1 = sin(2*pi*t);
% Calculate the norm of s1
s1_norm = norm(s1);
display(sprintf('L2 norm of s1: %.2f', s1_norm))
%% Create signal s2 as the FFT of normalised s1
% Scale FFT by sqrt(N) - this is convention used in MATLAB
s2 = fft(s1)/sqrt(length(t));
% Calculate the norm of s2
s2_norm = norm(s2);
display(sprintf('L2 norm of s2: %.2f', s2_norm))


What to do for the right FFT in Matlab (two peaks and incorrect amplitude)?

The somehow symeetric aspect of the complex DFT might be misleading, so in the past I shared a simple Matlab function FFTR.m for real signals only, that might yield the result you'd expect ($2$ magnitude at peak frequency):

as long as you don't try bordeline endeavours. That requires a little normalization and a special care at DC and Nyquist. The code could be:

% Signal definition
Amp = 2;
freqHz = 10000;
fsHz = freqHz*2+1;
dt = 1/fsHz;
sine = Amp*sin(2*pi*freqHz*(0:dt:1-dt));

% Graphic part
[fftR,fftAxe] = FFTR(sine',dt);
plot(fftAxe,fftR,'o-');
xlabel('Frequency')
ylabel('Magnitude')


Interpolation of magnitude of discrete Fourier transform (DFT)

First a demonstration that the squares of both
$$\begin{align}&[\dots, 0, 0, 1,\hphantom{-}1, 0, 0, \dots] \text{ and}\\
&[\dots, 0, 0, 1, -1, 0, 0, \dots]\end{align}$$

equal

$$[\dots, 0, 0, 1, \hphantom{-}1, 0, 0, \dots]\hphantom{\text{ and}}$$

but the squares of their sinc interpolations differ (Fig. 1):

Figure 1. Squares of sinc interpolations of $[1, 1]$ (blue) and $[1, -1]$ (red).

This demonstrates that it is not in general possible to recover the square of a band-limited signal from its uniform samples taken at the critical sampling frequency of the band-limited signal.

Let's test different interpolation approaches in Octave. The gold standard of band-limited interpolation is DFT–zero-pad–DFT:

>> format free
>> x = [1 2 3];
>> y = [1 2 3  0 0 0];
>> abs(fft(x))
ans =

 6 1.73205 1.73205

>> abs(fft(y))
ans =

 6 4.3589 1.73205 2 1.73205 4.3589

The last set of numbers are the magnitude values calculated from perfectly interpolated frequency domain bins. Let's try interpolating magnitude instead:

>> fft(fftshift(horzcat([0 0], fftshift(ifft(abs(fft(x)))), [0])))
ans =

 6 4.57735 1.73205 0.309401 1.73205 4.57735

Seems quite a bit off. Now let's try interpolating squared magnitude:

>> fft(fftshift(horzcat([0 0], fftshift(ifft(abs(fft(x)).^2)), [0])))
ans =

 36 25 3 -8 3 25

>> sqrt(ans)
ans =

 (6,0) (5,0) (1.73205,0) (0,2.82843) (1.73205,0) (5,0)

It's not only off, but also sqrt() returned a complex number for a negative interpolated value. So is the only valid way to interpolate the bin values?

Let's give this one more chance, by trying to interpolate frequency domain data upsampled by a factor of 2. Because of the even length transform, this requires duplicating the "Nyquist sample", so apologies if the code is getting hard to read.

>> z = [1 2 3  0 0 0  0 0 0  0 0 0];
>> abs(fft(z))
ans =

 6 5.55485 4.3589 2.82843 1.73205 1.77302 2 1.77302 1.73205 2.82843 4.3589 5.55485

The above is what we want. Let's try interpolating 2x upsampled magnitude:

>> fftshift(horzcat([0 0 0], fftshift(ifft(abs(fft(y)))), [0 0 0]))
ans =

 3.36365 1.10447 0.318175 0 0 0 0 0 0 -0.208949 0.318175 1.10447

>> ans(4) = ans(length(ans)-2)
ans =

 3.36365 1.10447 0.318175 -0.208949 0 0 0 0 0 -0.208949 0.318175 1.10447

>> fft(ans)
ans =

 5.79105 5.59483 4.56785 2.7273 1.5231 1.76882
 2.20895 1.76882 1.5231 2.7273 4.56785 5.59483

That is still off. Let's try interpolating 2x upsampled squared magnitude:

>> fftshift(horzcat([0 0 0], fftshift(ifft(abs(fft(y)).^2)), [0 0 0]))
ans =

 14 8 3 0 0 0 0 0 0 1.18424e-15 3 8

>> ans(4) = ans(length(ans)-2)
ans =

 14 8 3 1.18424e-15 0 0 0 0 0 1.18424e-15 3 8

>> sqrt(fft(ans))
ans =

 6 5.55485 4.3589 2.82843 1.73205 1.77302 2 1.77302 1.73205 2.82843 4.3589 5.55485

Now it works perfect! The take home message is to upsample (time domain zero pad) at least by a factor of two before trying to interpolate in the frequency domain, and to interpolate squared magnitude rather than magnitude. It works because taking the squared magnitude is the same as multiplying each bin value by its complex conjugate. Complex conjugation preserves the bandwidth of the bandlimited function represented by the data, so the multiplication doubles the "time domain bandwidth" because it is equivalent to time domain convolution. Note that when choosing the interpolation method the 2x upsampled squared magnitude is still critically sampled, so further oversampling should make accurate interpolation much easier.

I about answered my own question, but will accept further insight as an answer!

P.S. Just found out there is also interpft, which does the interpolation with less syntax.

Utilizing additional information

Interpolation becomes easier or even exact with additional information about the data, for example that it is a critically sampled squared time-shifted sinc. In that case, given the two samples $\alpha$ just before and $\gamma$ just after the largest sample, the time $-1<d<1$ of the peak can be calculated by:

$$d = \begin{cases}0&\text{if }\alpha = \gamma,\\
\frac{1 - \sqrt{1 - p^2}}{p},\, p = \frac{γ - α}{α + γ}&\text{otherwise,}\end{cases}$$

with $d = 0$ meaning that the sinc is shifted to exactly the time of the largest sample. The same interpolation can be done at half sample rate, which is the critical sampling frequency of the underlying function the magnitude of which equals that of a time-shifted sinc, with the two successive largest-valued samples $\alpha$ and $\beta$ of the square of the absolute value of the underlying function:

$$0 \le d \le 1\\
d = \begin{cases}0&\text{if }\beta = 0,\\
1&\text{if }\alpha = 0,\\
\frac{1 + p - \sqrt{1 - p^2}}{2p},\, p = \frac{\beta - α}{α + \beta}&\text{otherwise,}\end{cases}$$

The formulas are not affected by amplitude scaling of the data. For frequency estimation you'd rarely have a purely real time-shifted sinc, but if you do, there are exact interpolation formulas for it.


In pulse radars you perform matched filtering in receiver where you compute cross-correlation of the received waveform and reference one. Due to Doppler shift matched filter performance  will decrease. It's like decreasing of correlation function with increasing of time delay between waveforms. Doppler tolerance defines the maximum available Doppler shift for given waveform such that you can still achieve correlation peak bigger than some threshold (SNR is assumed high enough to produce such a correlation).

There is fundamental property of all waveforms used in Radar (and other apps): ambiguity function. Look https://en.wikipedia.org/wiki/Ambiguity_function. In a few words this function determines how correlation between 2 identical waveforms changes with time delay and also with doppler shift. Some waveforms have very narrow ambiguity function, while others - more spreading.

For example maximum length sequences have narrow AF both in time and frequency shifts. LFM has slightly wide AF especially in frequency shift. They are actually very tolerable to Doppler shift. COZAC sequences (used in LTE) have narrow AF in time axe but wide in frequency one.

So if you are choosing waveform for your application you should refer to its AF at first.


Frequency and phase range of the FFT of an image

The FFT calculates the discrete Fourier transform (DFT) of an image.
Since it is a separable transform (calculated on rows and columns independently) the explanation is the same for 1D signals.

The DFT of a signal, $X[k]$, is the projection of the signal $x[n]$ onto a complex frequency: $X[k]=\sum_n x[n]e^{-2\pi jnk/N}$.

$k/N$ can be considered to be the frequency. Since $k=0...N-1$ it starts with the lowest frequency (DC or the mean of the signal, $e^0$) and gets the highest frequency - a basis function $[-1,1,-1,1,-1,1,...]$, on the central frequency, $k=N/2$: $e^{2\pi j \frac{1}{2}n }=e^{j\pi n}$. In your example the result has been shifted (using fftshift) so that the lowest frequency is in the middle of the image, and the highest frequencies are in the border.

So, the range is the spatial frequency $0..N-1$ for each dimension. Since the result is complex you can present it using a magnitude image and a phase image, but the axes are the same. A certain point in the 2D transform image $(u,v)$ represents a frequency $e^{2\pi j(xu/N+vy/M)}=e^{2\pi j ux/N}e^{2\pi j vy/M}$ which turns out to be a spatial frequency with a certain direction and "speed".

E.g., $X[u,v]$ where $u=5, v=10, N=M=80$ represents the projection of the image onto the following basis function (here shown is the real part):

In your image, you can see most of the frequencies are concentrated near the origin (i.e. low frequencies), because most of the image is flat. Since there are edges (high frequencies) you can see some high magnitude lines on the transform image.
You can compare it to an image with high frequency content, such as white noise, in which there will be much more high frequency content.


How to read off power at different frequencies

I know my x-axis gives frequency in Hz.

That is not the case.

A DFT is a mathematical transform between discrete signals in time and frequency domain. Just as as the number of your sample doesn't give you time in seconds, the number of the DFT bin doesn't give you frequency in Hertz.

Instead, the number of the sample gives you time in "sampling intervals $T$", and thus, the bin number in the $N$-point DFT gives you frequency in "Nyquist rate fractions" $\frac {T^-1}N = \frac{f_{sample}}N$, whereby you need to pay attention – the DFT should be understood as something that shows you both negative and positive frequencies, so therange is $-\frac{f_{sample}}2$, to  $+\frac{f_{sample}}2$ (for odd and even lengths, take care of the boundaries when calculating what the frequencies of the DFT bases are). The zeroth DFT base is always the average/sum/DC of the signal.

  if possible in dBm.

Not possible, unless you know the proportionality between your digital input signal to physical power, ie. you've calibrated your measurement setup and know which power a digital $1+0j$ has. Then: Parseval's theorem is your friend.

Note that the DFT is nothing but a linear operator – a matrix that maps $\mathbb C^N\mapsto \mathbb C^N$; the question how to interpret complex frequency domain values is absolutely up to your understanding of your system. Sometimes, the values' magnitude squares can be interpreted as estimate for the Power Spectral Density at the frequency of the bin – at other times, that is not the case (especially if you don't know whether your signal is weak-sense stationary).

Sometimes, the bin of a DFT of a system is an impulse (ask a physicist about the meaning of the Fourier transform for solid state electronics), sometimes it's something like a pseudo-angle (ask an antenna designer), sometimes, it's proportional to speed of something (ask a radar engineer), sometimes, it's really just that — a complex, $N$-dimensional vector, that has no physical significance.

I'd like to quoute:

  One can Fourier transform anything—often meaningfully.
                                                          ~John Tukey


Image Matrix Vector Representation for the Degradation Model

Indeed since both expressions are scalars then they are equal to each other since the transpose of a scalar is the same scalar.
See in MATLAB as an example (Calculating $ {x}^{T} H y $ and $ {y}^{T} {H}^{T} x $:
>> vX = randn(10, 1);
>> vY = randn(10, 1);
>> mH = randn(10, 10);
>> vX.' * mH.' * vY

ans = -0.8618

>> vY.' * mH * vX

ans = -0.8618

As you can see, indeed both expressions are scalars.
As expected $ {x}^{T} H y = {\left( {x}^{T} H y \right)}^{T} = {y}^{T} {H}^{T} x $.


Compute a DFT of your signal over a period of $T$ samples, extract the phase $\phi$ (in the $[0, 2\pi[$ range) of the first harmonic, and use the point $\frac{2\pi-\phi}{2 \pi}T$ as your landmark.

This will not correspond to a particular characteristic landmark in the time-domain waveform; but the procedure to retrieve it is by design invariant to phase shifts, distortion adding higher harmonics, and to a large extent noise.

If the waveform of your signal is known in advance you can also use it as a template, and select as a landmark the lag, in samples, at which the cross-correlation between your observed period and your template is maximal.

Python code here:

import numpy
import pylab

def shift(s, k):
  n = s.shape[0]
  return s[(numpy.arange(n) + k) % n]

T = 256
t = numpy.arange(0.0, T) / T

original = numpy.sin(t * 2 * numpy.pi) + numpy.sin(3 * t * 2 * numpy.pi) / 3
time_shifted = shift(original, 40)
distorted = shift(original + original ** 3, 200)
noisy = shift(original + numpy.random.randn(1, T).ravel() * 0.2, 180)

signals = [original, time_shifted, distorted, noisy]

for i, s in enumerate(signals):
  pylab.subplot(411 + i)
  pylab.plot(s)

  # Phase of first harmonic method.
  first_harmonic_coefficient = numpy.fft.rfft(s)[1]
  angle_2pi = -numpy.angle(first_harmonic_coefficient) + 2 * numpy.pi
  angle_2pi = angle_2pi % (2 * numpy.pi)
  landmark = numpy.round(angle_2pi / (2 * numpy.pi) * T)
  pylab.plot(landmark, 0, 'r+')

  # Cross-correlation method.
  xcorr = numpy.fft.irfft(numpy.fft.rfft(s) * numpy.fft.rfft(original))
  landmark = numpy.argmax(xcorr)
  pylab.plot(landmark, 0, 'g*')

pylab.savefig('plot.png')


How to use MATLAB's $\tt xcorr$ to find index of a section inside a whole recording?

Tried with taking absolute value of the maximum ? I have updated the code, see if that works for you...

[ff,ind]=max(abs(c)); % finding the index of the offset

M=A(ind:ind+(length(B))-1); % the selection from A based on the new index

plot(B);
hold;
plot(M);


Kalman Filter - Updating the Covariance Matrix Step

Given the Covariance Matrix $ {P}_{k - 1 \mid k - 1} $ then:

$$ {P}_{k \mid k} = {F}_{k} {P}_{k - 1 \mid k - 1} {F}_{k}^{T} + {Q}_{k} $$

Where $ {F}_{k} $ is the Model Matrix at iteration $ k $ and $ {Q}_{k} $ is the Process Noise Covariance Matrix at iteration $ k $.


You can equalize magnitude and phase simultaneously by defining a desired complex frequency response
$$D(\omega)=M(\omega)e^{j\phi(\omega)}\tag{1}$$
with magnitude $M(\omega)$ and phase $\phi(\omega)$ chosen such that they compensate for the given magnitude and phase distortions.
An FIR filter approximating $(1)$ can be designed by using the following error function:
$$E=\sum_kW_k\big|D(\omega_k)-H(\omega_k)\big|^2\tag{2}$$
where $\omega_k$ are frequencies on a dense grid, $W_k$ are some non-negative weights allowing to emphasize certain frequency regions compared to others, and $H(\omega)$ is the frequency response of the FIR filter to be designed:
$$H(\omega)=\sum_{n=0}^{N-1}h[n]e^{-jn\omega}\tag{3}$$
Minimizing $(2)$ with respect to the filter coefficients $h[n]$ can be achieved by solving a system of linear equations. I've written a Matlab/Octave function that does just that: lslevin.m
Note that it may be necessary to add a linear phase (i.e., a delay) to the required $\phi(\omega)$ in order for a causal filter to better approximate the desired complex frequency response. As a guideline, the average delay implied by the desired phase should be about half the chosen filter length $N$. It is worthwhile experimenting with different additional delays because they may have a significant influence on the resulting approximation error.


What is so special about this filter coef = [1,2,1]?

It is a symmetric odd-sized FIR smoothing kernel, belonging to the class of Pascal or binomial filters that somehow sample a Gaussian kernel. Plus, its coefficients are simple dyadic integers, that can be implemented as bit-shifts 1/4 1/2 1/4. The coefficients sum to one, hence it is unit gain at DC.

In simpler word: (one of) the simplest real smoother preserving symmetric feature location. It's efficiency is however limited


Decimation of signal. How to plot frequency spectrum

The root of your problem is that you are not familiar with the Nyquist Sampling Theorem.

First of all, I question your premise that your sample rate is 70 MHz and that you can represent your signals, which range from 0 to 70 MHz, at this rate.  The Nyquist sampling theorem indicates that you can only represent signals up to half the sample rate, which is 35 MHz in your case.

When you downsample by a factor of four you are reducing your sample rate by a factor of four, which in turn reduces the frequencies that you can represent by a factor of four.  You are thus aliasing your signals, causing them to all be at the same frequency.  Thus, the one signal that you see after you downsample is really all four signals packed into one slot.


Notation confusion -- What is the correct operator for computation of the log-likelihood expression for complex valued data?

Without seeing your overall signal model we cannot really help. I suppose $y$ is your estimate? And hence, $L=L(y)$. Where is $h$ in your Likelihood equation? However, some general observations:

As Batman pointed out: likelihood corresponds to probability, which is necessarily a real value. So, if your equations yield complex values, you did something wrong. (Note that, depending on numerical accuracies etc. it can happen, that there is a very small imaginary part ($\approx 10^{-16}$) in your numeric computations. This part can be ignored.)
Commonly (but not always!) when in the real-domain you have a transpose in your equation, it changes to conjugate transpose (i.e. Hermitian operation) in the complex domain.
Noise variances commonly remain real-valued. However, noise realizations will be complex-valued.
You should understand the expression in your exponential as the squared distance between the measurement and the estimate. In real domain, the squared distance of $x,y$ is equal to $(x-y)^2$. However, in the complex domain, it becomes $\|x-y\|^2=(x-y)(x-y)^*$.


Which unit does the power spectrum of microphone output have?

This depends a little on how thorough you want to be. Your microphone measures pressure and not actual power. That means your time domain signal has the units of Pascal and the spectrum has Pascal squared per Hz. In order to turn that into Intensity you would have to normalize by the free field impedance of air which is density*speed_of_sound and roughly 410 Newton seconds per cubic meter (at sea level & room temperature). Then the units would be Watt per square meter per Hz. You would only get to Watts per Hz if you integrate over a full sphere . This requires either an omnidirectional source or a rather tedious measurement involving a large number of microphone positions.

There is added complexity if you want to calibrate in dBSPL (sound pressure level). The 0 dB reference for this is 20 micro Pascal (which is close to the human threshold of hearing) and a pressure of 1 Pascal corresponds to roughly 94 dBSPL.


The first peak is the DC component of your spectrum which is large compared to the AC components (your signal doesn't have negative value and its always above zero which lead to large DC bias)

Find average of your signal over time and subtract it from your signal to remove the DC component.

Good luck.


Higher side lobe levels in linear array beam pattern

There are several practical items that may cause reduced beamforming performance and lead the higher sidelobes. You haven't said how high exactly these sidelobes are, nor what levels you expected.
Some common areas to investigate further:

Array element spacing. To prevent the presence of grating lobes your element spacing needs to be less than $\frac{\lambda}{2}$. You haven't indicated the bandwidth of your signal, so I'll use the
400 KHz you gave and assume a speed of sound to be 1500 m/s. Then your spacing needs to be less than $\frac{1500}{2\cdot 400\cdot10^3}=1.875$ mm. Of course you should use the highest frequency of signal of interest for this calculation.
I'm assuming that you are using IQ sampling. Phase errors in the IQ demodulation signal or slew in the clocking (synchronization) across the 80 channels may degrade the performance. At this high frequency you will need very small tolerances to limit the effect of these phase errors.
Have you considered the effect of Multipath propagation? The effect of the signals bouncing off the side of the tank or the water surface can lead to signal levels with a higher intensity (constructive interference) than the level of your direct path.


Yes.  A high-pass filter removes low frequencies, and a constant baseline is the lowest frequency possible ("DC" = 0 Hz), so it will be completely removed and the baseline will become 0.


How do I manually plot the frequency response of a bandpass Butterworth filter in MATLAB without freqz function?

We know that in general transfer function of a filter is given by:

$$H(z)=\dfrac{\sum_{k=0}^{M}b_kz^{-k}}{\sum_{k=0}^{N}a_kz^{-k}} $$

Now substitute $z=e^{j\omega}$ to evaluate the transfer function on the unit circle:

$$H(e^{j\omega})=\dfrac{\sum_{k=0}^{M}b_ke^{-j\omega k}}{\sum_{k=0}^{N}a_ke^{-j\omega k}} $$

Thus this becomes only a problem of polynomial evaluation at a given $\omega$. Here are the steps:

Create a vector of angular frequencies $\omega = [0, \ldots,\pi]$ for the first half of spectrum (no need to go up to $2\pi$) and save it in w.
Pre-compute exponent $e^{-j\omega}$ at all of them and store it in variable ze.
Use the polyval function to calculate the values of numerator and denominator by calling: polyval(b, ze), divide them and store in H. Because we are interested in amplitude, then take the absolute value of the result.
Convert to dB scale by using: $H_{dB}=20\log_{10} H $ - in this case $1$ is the reference value.

Putting all of that in code:

%% Filter definition
a = [1 -0.5 -0.25]; % Some filter with lot's of static gain
b = [1 3 2];

%% My freqz calculation
N = 1024; % Number of points to evaluate at
upp = pi; % Evaluate only up to fs/2
% Create the vector of angular frequencies at one more point.
% After that remove the last element (Nyquist frequency)
w = linspace(0, pi, N+1);
w(end) = [];
ze = exp(-1j*w); % Pre-compute exponent
H = polyval(b, ze)./polyval(a, ze); % Evaluate transfer function and take the amplitude
Ha = abs(H);
Hdb  = 20*log10(Ha); % Convert to dB scale
wn   = w/pi;
% Plot and set axis limits
xlim = ([0 1]);
plot(wn, Hdb)
grid on

%% MATLAB freqz
figure
freqz(b,a)

Original output of freqz:

And the output of my script:

And quick comparison in linear scale - looks great!

[h_f, w_f] = freqz(b,a);
figure
xlim = ([0 1]);
plot(w, Ha) % mine
grid on
hold on
plot(w_f, abs(h_f), '--r') % MATLAB
legend({'my freqz','MATLAB freqz'})

Now you can rewrite it into some function and add few conditions to make it more useful.

Another way (previously proposed is more reliable) would be to use the fundamental property, that frequency response of a filter is a Fourier Transform of its impulse response:

$$H(\omega)=\mathcal{F}\{h(t)\} $$

Therefore you must feed into your system $\delta(t)$ signal, calculate the response of your filter and take the FFT of it:

d = [zeros(1,length(w_f)) 1 zeros(1,length(w_f)-1)];
h = filter(b, a, d);
HH = abs(fft(h));
HH = HH(1:length(w_f));

In comparison this will produce the following:


I see you've already accepted one of the answers, but you (or others) might be interested in: R. J. Webster, “Leakage Regulation in the Discrete Fourier Transform Spectrum,” Proc. IEEE, vol. 68, no. 10, Oct. 1980, pp. 1339-1341. -- Among the statements:

"A point insufficiently emphasized in the literature is that leakage in the DFT spectrum is a phenomenon of calculation.  The continuous distribution of leakage energy is determined only by the form of the window, and the DFT is a kind of lattice through which this continuous function is observed.  Because the lattice is fixed by convention, and the leakage distribution shifts according to the frequency of the signal, the discrete leakage components become dependent on the signal frequency."

(regarding the familiar 'picket-fence' effect): "The reader should note, importantly, that this effect is not a property of the harmonic nature of (the signal frequency fzero), as seems often suggested, e.g., [5], [6], but of the essentially arbitrary way G(f - (the signal frequency fzero)) is evaluated.  The apparently special relationship between harmonic frequencies and the rectangular window, in which the zeros of the window transform coincide precisely with the DFT sampling frequencies, is neither significant nor unique."

The author then goes on to show an example of using a shifted sum-of-cosines window centered on f = 10.25 that reduces the leakage to zero beyond bin 11 of his FFT when given an input signal of that frequency.


Purpose of image feature detection and matching

Image keypoints are a key feature in many Image and Video processing softwares, both industrial and academic. The principle behind is always the same:

detect some meaningful points in some images;
[optional] compute a stable description of the image part surrounding each keypoint;
match keypoints from an image (the template) to another (the query).

Now, some more detail.

Why detecting points? Basically, because a point is an intersection of two lines. As such, it is detected much more accurately (and in a more stable way) as lines or other features (area, etc.). This has been common knowledge in Computer Vision later proved in a paper by Mikolajczyk. This is rooted deeply in any point detector: for example, Harris point detector relies on the Hessian matrix to check that there are two crossing lines.

Why is the description part optional? Basic point detectors (Harris, SUSAN, Moravec, FAST) do not come with a descriptor. The process of matching them is purely location-based. Thus, they are mostly used for video processing. The intuition behind is that you can have many many points, then a robust estimation process (e.g., RANSAC). Recent keypoints however (SIFT and the ones that came after) embed an image patch descriptor, because they are meant to be used in more complex situations (object detection, large baseline matching...).

What's he difference between Harris/Hessian and SIFT/SURF? Besides the descriptor part, SIFT/SURF also include a localization in scale. Thus, when you detectd a SIFT/SURF point, you expect to find it again even if your image has been zoomed in/out, etc.

What's the difference between SIF and SURF? Well, SURF is actually mostly an approximated version of SIFT that is better suited for fast computations (by using integral images). Consequently, SURF is not as stable (as "good") as SIFT under affine transforms, but it comes at something like 1/10th of the computational cost.

Are they suited for realtime software? Yes, fore sure for Harris. SURF was designed to be fast, so a good implementation should be fast too. SIFT was not designed for speed, but you can find GPU implementations. But if you really want to be fast even on lower-end devices, then you need to check recent points such as FAST (detector), BRIEF, ORB, BRISK, FREAK (binary descriptors).

What can you do with them? Well, many things. By tracking points in videos then you can stablize them, augment them with objects, track objects. By matching points between cameras, you can automate the calibration process of 3D reconstruction pipelines. By matching points on parts of object, you can detect objects. And probably much more...


Is there a relationship between Average energy of transmission and power saving

First of all, 8 b/s/Hz is not a bit rate, but rather a spectral efficiency.

Second, you don't transmit $\log_2(M)$ bits per second, but rather $\log_2(M)$ bits per channel use, which is roughly $T_s$ seconds, where $T_s$ is the symbol time. So, the bit rate is $\frac{\log_2(M)}{T_s}$.

If you adjust your wording to the above by saying 8 bits/channel use, then, yes, you are right, you save energy by using MIMO, but the saving is $170/20$ not the reciprocal.


Sample & Hold: Estimate jitter delay of an ADC

I think the jitter to SNR formula is based on a gaussian jitter.

You use "rand" which yields uniformly distributed numbers. You should use randn() instead which yields numbers distributed according to a gaussian distribution.


A low pass filter has a frequency response that satisfies

$$|H(\omega)|\approx 0,\quad |\omega|>\omega_c\tag{1}$$

where $\omega_c$ is the cut-off frequency. A complex low pass filter must also satisfy

$$ H(\omega)\neq H^*(-\omega)\tag{2}$$

which causes its impulse response to be complex-valued. So the frequency response of a complex low pass filter is characterized by Eqs $(1)$ and $(2)$.

In discrete-time, you have the usual convolution sum

$$y[n]=\sum_{k=-\infty}^{\infty}h[k]x[n-k]\tag{3}$$

where $y[n]$ is the output sequence, $x[n]$ is the input sequence, and $h[n]$ is the filter's impulse response. Note that all sequences in $(3)$ are generally complex-valued. If you split up (3) into its real and imaginary part you get

$$y_R[n]=\sum_{k=-\infty}^{\infty}\left(h_R[k]x_R[n-k]-h_I[k]x_I[n-k]\right)\\
y_I[n]=\sum_{k=-\infty}^{\infty}\left(h_R[k]x_I[n-k]+h_I[k]x_R[n-k]\right)\tag{4}$$

where the subscripts $_R$ and $_I$ denote the real and imaginary parts, respectively. From $(4)$ you see that you have to implement two filters, one with impulse response $h_R[n]$ and the other with impulse response $h_I[n]$, and that you have to filter the real part as well as the imaginary part of the input signal with both filters.

As for applications in digital communications, anywhere where you have in-phase and quadrature components, you can model these two signals as a complex signal, and any filtering on this complex signal requires a complex filter (or, at least, a filter with complex input and output). However, in practice you most often encounter two special cases:

The filter is actually real-valued but has a complex input and output, i.e. you simply filter the real and imaginary parts of the input signal independently with the same filter. In Eq. $(4)$ that means that $h_I[n]=0$. This is the case when you apply a low pass filter after complex demodulation.
The filter is complex-valued but the input signal is real-valued, so you filter one (real-valued) signal with two filters (the real and imaginary parts of the impulse response). In Eq. $(4)$ that means that $x_I[n]=0$. This occurs when you filter a real-valued passband signal with an analytic band pass filter before complex demodulation. (Note that here the filter is a band pass, not a low pass).

A complete complex filtering operation with complex-valued inputs and outputs, and a complex impulse response (as specified by Eq. $(4)$) is used for modelling a digital communication system using an equivalent discrete-time baseband channel. The corresponding baseband channel filter is generally complex-valued.


What is the transfer function of this block diagram

As pointed out in the other answers, your friend's approach is correct. However, his end result for $H(z)$ isn't.

With $q[n]$ the input signal of the delay element, you get $q[n-1]$ at its output, and the given difference equation for $q[n]$ follows easily. The output $y[n]$ can then be written in terms of $x[n]$ and $q[n-1]$. In the $\mathcal{Z}$-transform domain you get algebraic equations in terms of $X(z)$, $Q(z)$, and $Y(z)$, from which you can eliminate $Q(z)$, leaving you with an expression for $H(z)=Y(z)/X(z)$:

$$H(z)=\frac{Y(z)}{X(z)}=-\frac23+\frac{z^{-1}}{1+\frac12 z^{-1}}=-\frac23\frac{1-z^{-1}}{1+\frac12 z^{-1}}\tag{1}$$

So it was basically just a sign error that messed up the end result.


Compare the quality of IQ data from different Software Defined Radios (SDRs)

First and foremost, I would recommend against over the air testing for this given the significant challenge in really being able to provide the same signal to each radio (since you have both temporal and spatial constraints that you cannot simultaneously meet).

I would instead use one GNU radio as a transmitter (or any other repeatable high quality source) to provide a consistent and repeatable waveform that can be used for each of the radios under test with a cabled test where the attenuation can be carefully and repeatably controlled over a wide range.

With this, there are several tests that immediately come to mind:

An ideal test waveform would be band limited complex Gaussian noise (pseudo-random noise as a repeatable pattern) since it will best occupy all possible positions on the complex IQ plane as well as all frequencies in the band of interest. This would then test the radio for a wide range of modulations, but if you have specific modulations in mind then those such waveforms for test purposes would also be ideal. (OFDM waveforms are well represented as complex Gaussian distributed bandlimited waveforms by the way, but it would be sufficient enough to simply create a complex random waveform and band-limit it for this generic case.)

An ideal test metric would be the correlation coefficient, $\rho$, as it will show errors from all sources including noise figure, quadrature error, LO phase noise, ADC quantization noise, IP3, IP2, etc.  For more information on that see Noise detection and How can I find SNR, PEAQ, and ODG values by comparing two audios?

The test should be done over a wide range of input power levels to show the related sensitivity and maximum power handling for comparative receivers. Thus the test metric for all receivers would be $\rho$ vs $P_{in}$.

The $\rho$ may be limited by the transmitter if that used is not of sufficient quality, so care must be taken to confirm that the $\rho$ for the transmit waveform is better than that which would result for all the tests performed (or represent the ceiling of the measurement). For the same reason start with the highest test power with known confirmed $\rho$ and adjust the power level in the cabled test receiver using an inline passive variable attenuator (with a rated much higher input P1dB) since any active device or variable amplifier is more likely to introduce additional distortion.

Additional tests involve how the receiver rejects interference such as intentional jammers or co-site transmitters. This is usually done with one tone testing in an adjacent band to help distinguish those receivers with better LO phase noise, and two tone interference testing in such a way to test for IP3 (and IP2 for zero-IF receivers in particular). The same test as above with the correlation coefficient vs input power is repeated, but in the presence of the jammers which can be done with unmodulated tones. Different recievers can have different sensitivities based on their architecture (or mis-managed spurs from a poor design), so a robust test would sweep the test jammers over a wide frequency range. If the radios are intended for specific applications, then the interference tests can be limited to cases appropriate for their known environment.

Instantaneous dynamic range in band would also be of interest.  As a simple test for this an in band single tone could be used at maximum input and the noise floor and spurious free dynamic range could be assessed in an FFT of the captured IQ signal. (You could also do specific linearity and single-tone and two-tone in band testing if you wanted to isolate these characteristics from what was otherwise in the catch-all $rho$ test described above).