Are the Daubechies 4 wavelet function and scaling function dimensionless?

Generally yes, discrete wavelets are often implemented as dimensionless filters, like discrete averagers or differences.
Additionally, wavelets are non dependent on the variable. If this is pertinent to your application, it is possible to equip the wavelet with dimensions: morally, the scaling coefficients ought to be of the same dimension as the variable, yet the detail coefficients can be thought as homogeneous to derivatives (at level 1), second derivatives (at level 2) etc. Consequently, with the proper normalization, they could be expressed in $\textrm{m}/\textrm{s}$, $\textrm{m}/\textrm{s}^2$, etc. according to each level.
If I detail the computation for the analysis/synthesis, decomposed onto the primal wavelet space ($\phi$ and $\psi$) and reconstructed in the dual wavelet space ($\tilde{\phi}$ and $\tilde{\psi}$):
$$f(x)=\sum <f,\phi_{j,k}>\tilde{\phi}_{j,k}(x) + \sum <f,\psi_{n,k}>\tilde{\psi}_{n,k}(x)\,.$$
The dimensional analysis entails that products or convolutions of $\phi$ and $\tilde{\phi}$, and of $\psi$ and $\tilde{\psi}$ should be unit-less. In other words, I would choose them  with "inverse" dimensions, and this could be done per scale.


Conditions for precoding matrix to preserve complex conjugate symmetry on DFT vector

I think that the entries in your matrix $\bf{T}$ must obey $a_{N-n+1,N-m+1} = a^*_{n,m}$.
This is saying that the entries in row $N-n+1$ are the same as the coefficients in row n but where the coefficients are conjugated and reversed. The pattern in $\bf{T}$ for $N=4$ is

$T_4 = \left[ \begin{smallmatrix} a_{11}&a_{12}&a_{13}&a_{14}\\ a_{21}&a_{22}&a_{23}&a_{24} \\ a^*_{24}&a^*_{23}&a^*_{22}&a^*_{21} \\ a^*_{14}&a^*_{13}&a^*_{12}&a^*_{11}\end{smallmatrix} \right]$

I'm sure someone will come up with a better and more precise answer.


For image quantization, is it necessary to sample image?

Quantization of a continuous signal is the process of representing an infinite set of input values, $x(t)$, that fall into an interval $I_k=[a_k,b_k]$, with a single level $y_k(t)$, the output, that represents all the inputs to that interval with some optimality criterion.

You can quantize an analog signal without sampling it. (however don't mix this with binarization of its amplitude). To do so, you can pass the analog signal from a set of reference thresholders, whose output will be constant as long as the input signal remains in a given interval. Binarization is not necessary if you do not want to convert this into a digital signal.


How transmission speed and bandwidth are linked?

Actually you also have to consider the noise floor of the communication channel and transmitted signal power to define the limit of data rate through an AWGN channel. The information capacity $C$ (bps) of the channel is given by the Shannon's formula:

$$ C = B \cdot \log_2( 1 + \frac{ \sigma_x^2}{ \sigma_n^2 } ) $$

where $B$ is the channel bandwidth in Hz., and $\frac{ \sigma_x^2}{ \sigma_n^2 }$ is the channel SNR (Signal to Noise Power Ratio).

One ideal case happens when there's no noise, then the capacity goes to infinity. This means that on a noiseless (ideal) channel you would transmit all the information (and more) instantly.


Why does an ideal high pass filter on this image generate gray values

Extreme positive values are black and extreme negative values are white. Values near zero are grey. Where there is not much change the high pass filter outputs values near zero. Hence the grey regions in the image.


Is the following system LTI (finite speed of signal propagation)?

No, if your system can (with sufficient accuracy) be modeled as a linear time-invariant system, then it is fully characterized by its impulse response. Since the frequency response is the Fourier transform of the impulse response, the system is also fully characterized by its frequency response. The frequency response does not only describe steady-state behavior but also transients (because it simply is a full description of the system). See also this answer to a related question.


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

It is not clear why you would want to do this, and perhaps you are trying to achieve something at a bigger picture level where this is the wrong approach. However to clear up your question:

Since your sampling frequency is 30MHz, your unique digital spectrum is only from 0 to 30 MHz (or -15MHz to +15MHz depending on which portion of the spectrum you want to take). The rest of the spectrum, out to 1880MHz and beyond is replicas of this spectrum stated above.

1880 MHz/30MHz is 62.6666... so the boundary at 30*62 = 1860 MHz is identical to what exists at 0, and then from there up to 1860 + 30 MHz = 1890 MHz is identical to what exists at 30 MHz. So I suspect from your approach above that the result was a frequency at 1880-1860 = 20 MHz, which is identical to the frequency at 1880 MHz!

This is because the following frequencies are all aliases of the same when the sampling rate is 30 MHz:

20 MHz, 20+30*1=50 MHz, 20+30*2 = 80 MHz, 20+30*3 = 110 MHz .... all the way to 1880 MHz, which is 20 + 30*62.

The following post on the unique span for any digital spectrum may help to further clarify this result:

https://stackoverflow.com/questions/44429592/what-is-the-result-inside-fft-bin-at-different-index/44434979#44434979

UPDATE: From the comment you said that you "want to see multiple signals at the same frequencies on the same spectrum". And I assume your spectrum is centered at 1880 MHz, and you may be simulating a multicarrier system in the 1850 MHz to 1910 MHz band (my guess) such as CA_2-4 or PCS. Assuming that is the case, then your maximum bandwidth of the signal you want to simulate is 60 MHz, and to simulate that I would recommend simulating using complex signals at
a sampling rate of 100 MHz (greater than 2x the single-sided BW of 30 MHz to appease Nyquist with additional room for evaluating effects beyond the band such as spectral regrowth or simulating filtering there, etc).  In this case 1880 MHz can just as easily be represented at DC, and another signal at 1181.25 MHz (for example) can be placed at 1.25 MHz in your simulation. Why? Because the items of concern for us with regards to the carrier is the complex envelope of the carrier only; with regards to anything in our system simulation, we don't really care about the actual high frequency cycle rate of the carrier itself- only the magnitude and phase of the envelope versus time.

Alternatively, if you want to simulate with a real carrier, this is easily achievable too by centering the real signal at $F_s/4$ (where $F_s$ is your sampling rate) in your digital frequency span. The frequency span needed to simulate a real 60 MHz BW is at least 120 MHz (leaving no room at band edges), and for that case I would likely use 200 MHz which would be more than sufficient. This would result in the "1880 MHz carrier" being placed at 100 MHz.

Here is a close example I have that demonstrates this which was for a multi-carrier CDMA system (reverse link).

This is actually the radio receiver design but I draw your attention to the 35 MHz Digital IF, that is then represented with a 20 MHz sampling rate. In the end this is providing a real digital output (not a simulation) of what is at the 35 MHz carrier, which is also the same signal at the 1880 MHz carrier. The digital output is sampled at 20 MHz and centered at 5 MHz, which is $F_s/4$. The reason to center a real signal at or close to $F_s/4$ is that it maximizes the distance in frequency to the images, which in practical implementation such as what is shown, it simplifies the analog filtering required. In a simulation it will center your signal in the unique digital span.

This drawing may also help illustrate the point of not needing to simulate every cycle of the actual carrier:

A bandwidth modulated signal is shown in time and in frequency. For the case of time we see it as a vector moving around on the complex I Q plane versus time, changing in magnitude and phase (which is what we are interested in!). It moves relatively slowly due to being bandwidth limited (as opposed to white noise which would instantaneously be at any other arbitrary point at every next instant in time; we would just see a 2 dimensional Gaussian cloud in that case). In frequency we see the signal is centered at 0, and has some bandwidth as given by the rate of modulation and any filtering or similarly pulse shaping that was applied.

This next slide demonstrating frequency translation also helps to make the point I am trying to make here. On the left we see the same signal except at a carrier frequency of $F_o$.  The IQ time domain plot in this case is showing how the phase and magnitude of the signal is changing versus time relative to the carrier.

Consider the carrier alone with no modulation: When the carrier is at 0 as in our first case, the IQ plot would be a vector of constant magnitude and angle, not moving versus time (If the angle was 0 and magnitude 1 then it would be e vector on the real or I-axis centered on the origin pointing to the right). When the carrier is a positive frequency 1880 MHz signal (and to simplify this description consider the signal to still be complex), the IQ diagram would now be a vector rotating counter-clockwise at the rate of 1880 MHz.  If we could observe this with a strobe-light flashing at the same 1880 MHz rate synchronized to when the vector passes $1\angle 0$, we would then see the same case as when the carrier was at 0: a vector centered on the origin pointing to the right and not rotating.  Consider the same thought process with the modulated signal as in the picture below. When the signal is at the 1880 MHz carrier, the IQ diagram would be spinning at the 1880 MHz rate, but if we could observe it with the strobe light as previously described, we would see how the magnitude and phase deviates from 1880 MHz over time at each sample of the strobe light, and it would look identical to the same signal centered at 0. Since the signal is changing so slowly compared to 1880 MHz, there is no need for the strobe light to be going that high of a rate, we could equivalently run the strobe light at some sub-multiple of 1880 MHz (such as 100 MHz) and still observe the same result with the signal centered at 0; as at the 1880 MHz rate the signal has changed insignificantly in magnitude and phase from the previous sample.


How to add white gaussian noise to an image in Matlab?

Additive noise is usually mean-free, i.e. the mean should be 0. A noise level of 10 dB = 1 B ("bel") usually means that the variance of the noise is by a factor 10¹ = 10 smaller than the variance of the image. As the documentation states, the variance parameter to imshow is interpreted under the assumption that the image data are doubles in the range [0, 1].

Here's an example:

Load sample grayscale image included with Matlab:

[I, map] = imread('eight.tif');

That image is of data type uint8, i.e. integers in the range 0 to 255. Let's convert it to a double in [0, 1]:

I = double(I) / 255;

Compute the image variance

v = var(I(:));

and tell imnoise to add Gaussian white noise with mean 0 and a tenth of that variance:

J = imnoise(I, 'gaussian', 0, v / 10);

This is the result:

Because imnoise clips the data to the range [0, 1] after adding noise, the actual noise variance is smaller,

var(J(:) - I(:)) / v

gives something like 0.093.


Derivation of Yule-Walker Equation for LPC

Yes, your derivation is correct.
There is just a minor inconsistency in your definition of the auto-correlation function $r(k)$. If you choose to define
$$r(k) = E[s(n-k)s^*(n)]$$
then $E[s(n-k)s^*(n-i)]$ equals $r(k-i)$ or, equivalently, $r^*(i-k)$, but not $r(i-k)$. So the last equation should read
$$\sum_{i=1}^pa^*_i r^*(i-k) = -r(k)$$
You take the derivative w.r.t. to $a_k$ because these are the coefficients that you can choose as to minimize the error. They are the unknowns in your problem.
The fact that $e^*(n)$ can be viewed as a constant when taking the derivative w.r.t. $a_k$ has to do with how computing the gradient w.r.t. complex variables works. It is equivalent to taking the derivatives w.r.t. to the real and imaginary part of the coefficients separately, but the complex technique is more elegant and easier (at least as soon as you understand it). I've explained this technique in this answer.
The complex conjugates just pop up because you consider $|e(n)|^2$ in your optimization objective: $|e(n)|^2=e(n)e^*(n)$. You wouldn't need them at all if you only considered real-valued coefficients $a_k$ and real-valued signals $s(n)$.


Reducing the number of data points of an FFT

Do you want to know specific frequencies? Go for the Goertzel algorithm.

Do you need to know every 2nd frequency? This corresponds to downsampling in the frequency domain, which is aliasing in time domain. You can do it like this (where I have reduced the sampling rate and signal duration for better plotting, obviously it also works with 30s and 500Hz):

Fs = 50 # 50 Hz for less data
x = np.random.randn(2*Fs)  # The audio signal: 2s for less data for better plotting here
L = len(x)

# Calculate full FFT for reference
X = np.fft.fft(x)
f1 = np.linspace(0, Fs, L, endpoint=False)

# Calculate every 2nd sample of FFT
x2 = x[:L/2] + x[L/2:] # Perform the aliasing operation in time domain
X2 = np.fft.fft(x2)
f2 = np.linspace(0, Fs, L/2, endpoint=False)

plt.plot(f2, abs(X2), 'go-')
plt.plot(f1, abs(X), 'rx-')

plt.xlim((0, Fs/2))

Obviously, the approach of aliasing also works for every Nth sample in general, if $L$ can be divided by $N$ without remainder.


Even i had confusion differentiating between these two terms at the beginning ,
have a look at this explanation from one of the Award winning DSP books.

Amplitudes,

Magnitudes,


Incorrect initial statement.  The sampling rate must be greater-than (not just at-least or equal-to) twice the highest frequency present in a band-limited base-band signal, for aliasing (and thus potential information mixing/corruption/ loss) not to occur.

The maximum frequency present is usually by design, by putting a (estimated  or assumed to be sufficient) low-pass filter before the sampling mechanism.  (Or the mechanism of the signal or sampler itself includes such a low-pass frequency response.)


Why the delta at the end of the approximation?

Here the signal is in continuous time domain. We can approximate any signal with weighted integral of unit impulse. As the signal is in continuous time domain integration is used instead of summation.

The actual equation is ,
$$x(t)=\int_{-\infty}^{\infty}{x(t_0)\delta(t-t_0)dt}$$

In your equation ${lim}_{\Delta\rightarrow0}\sum^{\infty}_{-\infty}$ stands for integration so $\Delta$ is required at the end, which stands for dt in integral.

Hope you clear that.


Proof that module of FT of 2 independent signals is sum of modules

Assume $x[n]$ and $y[n]$ represent two jointly WSS (wide-sense stationary) discrete-time random processes. Further assume that they have zero means, without losing the generality.

If $x[n]$ and $y[n]$ are independent, then they are also uncorrelated, and since they have zero means, they are also orthogonal; i.e.

$$ E\{ x[n]y^*[n-k]\} = 0 \tag{1}$$

Let $$z[n] = x[n] + y[n] \tag{2} $$ be the sum process, then PSD (power spectral density) of $z[n]$ is

$$ S_z(\omega) = FT\{ r_z[m] \} \tag{3} $$
where  $r_z[m]$ is the Auto-Correlation sequence of $z[n]$:

$$ r_z[m] = E\{ z[n]z^*[n-m] \} \tag{4} $$

Based on Eq.1 and Eq.2, expand Eq.4 as:

$$
\begin{align}
r_z[m] &= E\{ z[n]z^*[n-m] \} \\
& = E\{ (x[n]+y[n])(x^*[n-m]+y^*[n-m]) \} \\
& = E\{ x[n]x^*[n-m]\} + E\{ y[n]y^*[n-m]\} + \\
   &+  E\{ x[n]y^*[n-m]\} + E\{ y^*[n]x[n-m]\}  \\
& = r_x[m] + r_y[m] +  0 + 0 \tag{5}
\end{align}
$$

The zeros on the last line follow because of Eq.1, ($x[n]$ and $y[n]$ are othogonal).

Hence what follows is that the PSD of $z[n]$ is the sum of PSDs of $x[n]$ and $y[n]$ :

$$ S_z(\omega) = S_x(\omega) + S_y(\omega) \tag{6}$$

Everything up to this point is theory based on mathematical (probabilty) model of random processes. What you are after, however, is about the practical computation of the PSD estimation and how the independence of $x[n]$ and $y[n]$ is reflected into tho Fourier calculations.

Assuming we use the periodogram estimate for PSD :
$$ \hat{S}_x(\omega) = \frac{1}{N} |X(\omega)|^2 \tag{7}$$

then what we want to show is, as you have also worked out,

$$ \hat{S}_z(\omega) = \hat{S}_x(\omega) + \hat{S}_y(\omega) \tag{8}$$

which implies that

$$ |FT(x+y)|^2 = (X + Y)(X^*+Y^*) = |X|^2 + |Y|^2 + XY^* + X^*Y = |X|^2 + |Y|^2 \tag{9} $$

and Eq.9 implies that:

$$ XY^* + X^*Y = 0 \tag{9}$$

Eventually it's convincing to show that $XY^* = 0$.

To show this is you can use the ergodicty assumption. If $x[n]$ and $y[n]$ are ergodic in their cross-correlation, then the theoretical expectation can be replaced by the practical sums as:

$$E\{x[n]y^*[n-m] \} = \lim_{N \to \infty} \frac{1}{N} \sum_{n=0}^{N-1} x[n]^*y[n-m] \tag{10}
$$

then it can be shown that the estimate is equivalent to convolution of the sequences :

$$
\frac{1}{N} \sum_{n=0}^{N-1} x[n]y^*[n-m]  = x[n] \star y^*[-n]\tag{11}
$$
But from Eq.1 we know that the $x[n]$ and $y[n]$ are orthogonal implying that

$$ x[n] \star y^*[-n] = 0 \tag{12}$$

finally, if this sum is zero then it's Fourier transform  also is zero, completing the proof.

$$ x[n] \star y^*[n] = 0 \implies X[k]Y^*[k] = 0 $$

Note that, as MArcus has already commented, the equality holds only for the ideal case of infinite length process realization and for any practical finite length sample sequence, the Eq.s (6,7,8,9,10,12) won't hold true;i.e., they will be approximate).


How would I go about designing a difference equation for a given type of filter?

Matlab's filter() function implements a difference equation by calculating

$$\begin{align}a_0y[n]=&b_0x[n]+b_1x[n-1]+\ldots +b_Mx[n-M]\\&-a_1y[n-1]-a_2y[n-2]-\ldots-a_Ny[n-N]\end{align}\tag{1}$$

for every sample index $n$, where $y[n]$ is the output sequence, and $x[n]$ is the input sequence.
The coefficients $a_i$ and $b_i$ are stored in the vectors a and b. Note that normally the first element of the vector a is chosen as $a_0=1$. From (1) you see that in order to implement the difference equation you need to store past input and output values. This is usually done with a circular buffer.

Normally you don't talk about "designing a difference equation", but you design a filter, i.e. you compute the coefficients $a_i$ and $b_i$ in such a way that the filter performs some chosen task. There is a large variety of filter design algorithms available. Take a look at this answer for some basic background information on filter design (this is a huge topic in itself!).

In most cases you would use some filter design package (e.g. available in Matlab or Octave), design the filter (i.e. compute the filter coefficients), and then you can implement the filter according to Eq. (1). Writing your own filter design routine can be a nice exercise, but since filter design is a very mature topic, you'll basically be re-inventing the wheel. Implementing the filter is a different thing, and it can be useful to write your own implementation if you want to perform some real-time filtering.

Note that a direct implementation of Eq. (1) is only one of many ways to implement a filter. There are many possible implementations, all with their advantages and disadvantages. They all implement Eq. (1), but they perform computations in different orders, and they use different intermediate variables, which results in very different numerical behaviors. Of course, they would all be equivalent if it were possible to perform exact computations without any numerical errors. Especially for fixed-point implementations it is very important to choose an appropriate structure, because a direct implementation of Eq. (1) will be disastrous for larger filter orders (i.e. larger values of $M$ and $N$).

Note that Eq. (1) is a time-domain implementation. If you're interested in frequency domain implementations, search for "overlap-add" and "overlap-save" methods. However, I would start with implementing a time-domain filtering routine if you're working on your first implementation. Also note that for frequency domain filtering you need an FFT routine.

If you can use floating point arithmetic, I would as a first step try to implement Eq. (1) directly, and compare your results with Matlab's filter() function. For computing the coefficients you should use Matlab's filter design functions. If you then have any more specific questions you could come back and formulate a new question.


Yes, $a'(\phi)$ is a steering vector. In fact, the third equation (almost) describes a phased array that is capable of steering lobe patterns at different frequencies.

Each element of the array has a filter behind it which forms a weighted sum of (delayed) samples in Time. The same idea is applied one level above by now summing (delayed) samples of the received signal from the different elements that are in different points in Space.

In STAP, the adaptation in Time maintains the coherence of the probe signal (or carrier) against noise / doppler and other "structural" distortions and the adaptation in Space maintains the coherence of the signal along a specific direction by trying to minimise the covariance between different elements. Think of it as one convolution per element of the array(typical filtering in the Time domain) and one more ACROSS the same time instant (the output) of the element filters (Space).

BUT! With this configuration you can tune the Time filters at some sinusoid and then with the frequency of that sinusoid tune the Space filter to receive from a given bearing. To do the same thing simultaneously, you would need a Fourier Transform behind each element to decompose the received signal into a bank of sinusoids and then apply different space filters ($a'(\phi)$) to those outputs for specific bearings. In this way, you could have the array listening to some low frequency from the left and some high frequency from the right simultaneously. I feel that this is what is meant with "...the more generic case is that the array response on all channels is not equal (more like an array manifold) and is thus..."

Hope this helps.


Transform lowpass filter to highpass filter, what happened to unit sample response and difference equation?

I know that multiple a lowpass filter signal by (−1)n can transform it to a highpass filter.

Depends a bit what you mean by that. If you assume that this will turn 1 kHz lowpass signal into a 1kHz high pass signal, your assumption is plain wrong. What happens here is that SHIFT your signal frequency spectrum by half the sample rate, $f_s$. Due to the periodicity of the spectrum, this looks like the spectrum has been mirrored at $fs_/4$ (at least for real valued signals), i.e. DC become Nyquist, Nyquist becomes DC.
So this operation flips the entire signal spectrum, regardless of whether it's been high-pass, low-pass, pink, white, etc.
Loosely speaking, if you have signal that's been low passed and you flip it, the spectrum may look a bit like that of a high passed signal, but it's completely different operation than filtering the original signal with a complimentary high pass.

But exactly what happened to the unit sample response and difference equation? Assume that the original lowpass filter signal has a difference equation.

This is NOT an LTI operation, so difference equations and unit sample response do not apply here. Specifically, it's time variant.


Subtract background from 2D signal using spectral methods

There are certainly ways to "kill" the underlying signal but in this case, it is probably better to do it gracefully.

  I'm not even very clear about whether I should be looking to remove the highest frequencies, the lowest frequencies, or both.

In general, completely supressing a (usually) low frequency background oscillation is done by using a high pass filter.  Other ways this can be achieved is by applying a median filter to extract the variation and then subtract it from the original signal, or even peforming a low order polynomial fit to the data and then extract the surface this polynomial describes. For slightly more information on those, please see this response.

Here however, the background signal is not just any nuisance but it is related to the measurement. Therefore, if we just try to chuck it away, we may be introducing a problem to the measurements we are trying to get to. So, the better approach would be to see if it is possible to take the component out somewhat more reasonably.

More so, when you think about the unknowns, which are not many. In the "Beta" dimension, we have 1 "revolution". The same applies for the "Alpha" dimension, but there, we also have a "known" phase shift and so only observe part of the sinusoid. In the "Intensity" dimension, we have a mean background over which the counts vary. So, we know the frequency of the sinusoid as well as its phase (for the "Alpha" dimension) and we would be looking at its amplitude.

In 1 dimension, the "fit" looks like this (here in GNU Octave but easily translatable to other platforms too):

Fs=256; % Sampling frequency in Hz
T=1; % Duration of the signal in seconds
t=0:(1./Fs):(T-(1./Fs)); % Time vector
p=2.*pi.*t; % Phase vector

% Setting up the "signal"
average_count=30; % Self explanatory
variation_amp=6; % The amplitude of the sinusoidal variation
sig_amp=4; % The amplitude of what here simulates the signal.
variation_phase=pi./0.8; % The phase of the sinusoidal variation.

s = average_count + variation_amp.*sin(p+variation_phase) + sig_amp.*rand(size(p));

% At this point, we know nothing about the parameters, except
% the frequency of the sinusoid, because we know that due to
% the experimental conditions, there has to be sinusoidal
% variation proportional to an angle that we control.

signal_mean = sum(s)
signal_amp = sum(s.*exp(-j.*p));
v = (signal_mean + signal_amp.*exp(j.*p) + conj(signal_amp).*exp(j.*-p))./length(p);
% Recovery

u = s - real(v);

% Display

subplot(311);plot(s);grid on;title("Signal");xlabel("Discrete time (sample)");ylabel("Amplitude");
subplot(312);plot(real(v));grid on;title("Recovered variation");xlabel("Discrete time (sample)");ylabel("Amplitude");
subplot(313);plot(u);grid on;title("Adjusted signal");xlabel("Discrete time (sample)");ylabel("Amplitude");

This produces:

Where, basically, what we have recovered after subtracting the predictable part is that additive signal which here looks like noise but could be any remaining "form".

This is basically Discrete Fourier Transform (DFT) but just for 1 component, the one we know (or, can predict) is in the signal. These lines:

signal_mean = sum(s)
signal_amp = sum(s.*exp(-j.*p));

Perform the forward transform. And this line:

v = (signal_mean + signal_amp.*exp(j.*p) + conj(signal_amp).*exp(j.*-p))./length(p);

Performs the inverse transform but just for one component. Notice here that signal_amp recovers both the amplitude and the phase of the sinusoidal we are after.

A similar approach would be followed in two dimensions but with the necessary adjustments:

Fs=256; % Spatial sampling frequency in Lines per unit of length
T=1; % Length of image in units of length (assumed square here for simplicity)
t=0:(1./Fs):(T-(1./Fs)); % **Space** vector
p=2.*pi.*t; % Spatial phase vector
[Xp, Yp] = meshgrid(p); % Spatial phase vectors for each dimension.

% Setting up the "signal"
average_count=30; % Self explanatory
variation_amp=6; % The amplitude of the sinusoidal variation
sig_amp=4; % The amplitude of what here simulates the signal.
% The phases of the sinusoidal variation.
variation_phase_x=pi./1.2;
variation_phase_y=pi./0.8;

s = average_count + variation_amp.*(sin(Xp+variation_phase_x)+sin(Yp+variation_phase_y)) + sig_amp.*rand(size(Xp));

rec_average = sum(sum(s));
Ay = sum(s.*exp(-j.*Yp))
Ax = sum(s.*exp(-j.*Xp),2)

rec = rec_average./prod(size(s)) + (Ax.*exp(j.*Xp) + conj(Ax).*exp(j.*-Xp) + Ay.*exp(j.*Yp) + conj(Ay).*exp(j.*-Yp))./256;

subplot(311);title("Original");surf(s);xlabel("Dim X");ylabel("Dim Y");zlabel("Amp");
subplot(312);title("Component");surf(rec);xlabel("Dim X");ylabel("Dim Y");zlabel("Amp");
subplot(313);title("Adjusted");surf(s-rec);xlabel("Dim X");ylabel("Dim Y");zlabel("Amp");

This produces:

Where as you can see, the variation is gone and only the additive component has remained. The two snippets have an almost one-to-one correspondence to show both the similarities but also adjustments you have to apply when considering higher dimensions.

In your case, since you are dealing with "counts" that must remain positive, you can ommit the DC (the first component that is simply the average across the signal) and only remove the variation.

Hope this helps.


The autocorrelation function is the inverse Fourier Transform
of the power spectrum. I have no idea where you got your formula
from or what your calculation is
meant to be doing. Are you sure that you are reading the
book/journal article/application note/website/Wikipedia or whatever
correctly?

In any case, if you want to compute an autocorrelation
as an ifft, you should be careful about making sure that
the power spectrum vector is set up right. If the
autocorrelation function that you compute as the iffy
does not have a peak at the origin (or as the first
entry in the vector for MATLABi folks) and have the
symmetry property $R[i] = R[-i]$ in real life and
the appropriately translated version in MATLAB, something
is awry.


Comparison of Bilateral Filter and Anisotropic Diffusion

Your question is a bit harsh, because it's kind of vague. I will give you a few points, maybe it will help.

What's the same?

The intuitions behind both bilateral filtering and anisotropic diffusion are the same:

averaging is good to remove random noise;
averaging should only concern pixels that belong to the same region (in the sense that they are pixels of a given color of a given object):

to avoid artifacts at the borders between regions (most notably, blur)
to improve the denoising efficiency (without outliers inside, averaging is more efficient).

What's different?

Bilateral filtering and anisotropic diffusion will differ by how they achieve these two goals (under the constraint that an image segmentation is not available).

Usually, anisotropic diffusion is expressed in a variational framework, where some image functional is minimized. These functionals usually include a term that depends on the gradient of the image to be denoised, because this term allows to trigger the anisotropy property. The solutions to these problems are usually obtained iteratively via some gradient descent, but not always (in image morphology inspired frameworks such as mean curvature motion for example, image level lines are moved according to their curvature).
Basically, the solutions to these functionals correspond with isotropic diffusion inside homogeneous areas, and diffusion across image edges is stopped by a weight depending on the image gradient.

Bilateral filtering is a pixel-based approach. For a given pixel, its denoised counterpart is obtained by the weighted average of its neighbours, where the weights are given by some function that depends on their color similarity and image distance.
Unlike anisotropic diffusion then, the problem is not solved globally over the image but over each a neighbourhood surrounding each pixel.
Furthermore, there is no explicit gradient-based barrier at a given pixel: this effect will be more or less pronounced by the value of the decay put in the visual distance function.

Note that when the bilateral neighborhood of a pixel goes infinite, you have a new algorithm called nonlocal means.

And now.. the speed

In terms of quality, I would consider that either approach is good. Anisotropic diffusion has some caveat though:

since it's usually implemented as some gradient descent or PDE solver, a quick-and-dirty implementation can run into numerical issues;
in my opinion, it's easier to tune the decay of the distance function in bilateral filtering than the data term in anisotropic diffusion functionals.

Remains the speed... This part depends a lot on the practical implementation.

When bilateral neighborhood size gets large (OpenCV claims large is above 5 pixels) then bilateral filtering is slow.
You can use some tricks (Gaussian approximated by boxes, pre-selection criterion...) to accelerate the code. In fact, there's even a significant part of the literature on bilateral filtering that is dedicated to speeding it up.

For anisotropic diffusion, it depends in the numerical scheme implemented to solve the problem. Some are better (faster and more stable) than others: typically, TV regularized constraint is nowadays implemented with Chambolle's projectors (fast) or with Nesterov's acceleration (see FISTA -- faster).
Furthermore, these schemes are well suited to massively parallel implementations (GPGPU, data parallel multithreading) that reduce greatly the computation time.


Binary robust independent elementary features ("BRIEF") understanding

First of all, there's no such thing as 'template' in this paper - the word 'template(s)' has a different meaning in Computer Vision.

The method used in this paper is relatively straight-forward. Let me break it down for you. There are three important things that you need to do when doing tasks such as object recognition, image matching, image stitching, and similar others, using Local Features

The first important thing is Detection; in this step you want to detect points-of-interest or keypoints, and what that means is that you want to choose local points (basically small patches) that you think are interesting in the image, there are many ways to do that; this paper doesn't contribute in this area. However, it seems that they use SURF feature detector and CenSurE keypoints (you can look them up if you want to know more about them, I won't talk about this step except that they use features such as gradients and so, which means that if you have a photo, it's unlikely that a point in the middle of the sky would be chosen as an interest point, that's because the pixels around that point are all same intensities, things that are 'busy' tend to be chosen as interest points (e.g. table/building edges/corners)).

After detection is done, Feature Description follows. You know the interesting points in the image and now you want to describe them (basically you want to describe the points/patch around the interesting points). SIFT is one popular feature descriptor. In this paper, they come up with a new one, called BRIEF. BRIEF is based on comparisons, so let's say we have a patch (50 pixels by 50 pixels), we choose two points and compare the intensities of the two points, if the first point is larger than the second point, we assign the value '1', else '0', we do that for a number of pairs and we end up with a string of boolean values. Now the big question is 'how do you pick the pair of points?', in the paper they explained 5 ways, relatively similar, I will describe the first one. What you do is you uniformly (equal probability) choose a point from -S/2 to S/2, in our example we said that the patch size is 50, so we choose a point between -25 to 25. Assuming that the 0,0 coordinate is located at the center of the patch. So here's an example;

We want to select the first pair, each point consists of (X,Y) coordinates so we uniformly select the X-coordinate of the first point, then the Y-coordinate of the first point, let's assume it's (10,-1), now the second point; we uniformly select the X-coordinate for the second point, and the Y-coordinate for the second point, let's assume it's (-2,20), now we get the intensity values for each point and see which one has the larger intensity values - if first is larger we assign the first boolean value to be '1', if not, we assign '0'. We do that for many pairs, and we end up with a vector of boolean values.

*Very important point: I assume for this to work, you will have to specify a seed value before the random generator. What that means is that you want to choose the same values for each patch - this is super important because when you compare/match patches, the whole system will break down if they're not chosen the exact same way. - hopefully this will make sense when you read the matching step.

So, we do that for each interest point that was detected by the detector. For each interest point, we will get a vector of boolean values.

Now to match two images (third step; matching), we do the exact same thing for the other image, we detect, then describe using the BRIEF. For example let's say we have 10 interest points in each image (this can always work if we get the 10 most interesting points in each image), we use BRIEF to describe each patch using for example 50 pairs, so each image will be described by 10 vectors of 50 boolean values.

To compare the two images, we find the nearest neighbor of each vector from the first image to each other vector from the second image. We use Hamming distance which is pretty fast, example of hamming distance

hammingDistance((0, 1, 1), (0, 0, 0)) = 2

hammingDistance((0, 1, 1), (0, 1, 1)) = 0

hammingDistance((0, 1, 1), (1, 0, 0)) = 3

hammingDistance((0, 1, 1), (1, 1, 1)) = 1

Basically how many wrong correspondences

Hope that helps


For spatial filtering there are numerous beamforming techniques and it will be difficult to compare all here, there are numerous resource available on internet.

As related to your specified beam forming techniques,

Phased shift beamformer also known as conventional beamformer:
This beamforming techniques tries to add constructively the signal from the provide direction. This algorithm thus shifts the signal such that they align perfectly at the provided direction and thus sums up constructively. This algorithm has the worst angular resolution amongst different beam forming techniques as all signals from other directions are also included into the result. The resolution is in the range of the beam width of linear array.
Minimum variance distortionless response (MVDR) beamformer also known as CAPON beamformer: This beamformer is based on matching the array response matrix with the inverse of the sample correlation  matrix, thus enhancing the signal from considered direction and also suppressing the effects from the interfere signals from different direction other than the direction under consideration. The resolution of MVDR is far better than phased shift beamformer.

Among these two algorithm stated above the MVDR has benefit related to separation between two signals but on the expense of more computations.
It would be of worth to also look into the Linear constraint minimum variance beamformer (LCMV).

The implementation of MVDR and LCMV weights for beamformer is present in MATLAB MVDR weights and LCMV beamformer.


Dimension of channel matrix elements in MIMO systems

Sorry, I'm starting again.
$\mathbf{H(t,\tau)}=(h_{ij}(t,\tau))$.

If I have $\mathbf{H}(t,\tau)=\mathbf{H(\tau)}$ and $N_t=1$ then $\mathbf{H}(\tau) \in \mathbb{C}^{N_r \times 1}$


Gain function calculation (frequency response)

One way to see the equivalence is to first factor out the $0.5$ in your initial transfer function, i.e.

$$
H(f) = 0.5(1 + e^{-i2\pi f}).
$$

Then you can rewrite the $1$ and $e^{-i2\pi f}$ as follows:

$$
H(f) = 0.5(e^{i\pi f}e^{-i\pi f} + e^{-i\pi f}e^{-i\pi f}) = 0.5(e^{i\pi f} + e^{-i\pi f})e^{-i\pi f}.
$$

Your result then follows as a trivial transform.


Does Radio Communication Have to Account for the Doppler Effect?

Does radio communication have to account for the doppler effect?

Yes.
Would be pretty terrible if that wasn't the case: RADAR wouldn't work!

Do moving objects such as planes and rockets have to account for this,

Yes. Phones in cars and trains, too. Your 5G NR phone of the future operating above 60 GHz, will have to do that at walking speeds, too.

or is this movement negligible compared to the speed of light?

You can calculate that yourself! With $f_0$ being the original frequency of your signal, and $c_0=3\cdot10^8\,\frac{\text{m}}{\text{s}}$ the speed of light, the Doppler shift as seen by an object travelling at speed $v$ is $\Delta f= f_0\cdot\frac{v}{c_0}$.
So, no, for most microwave communication (that's basically everything "broadband" that you do these days – TV, cellular, Wifi, drone video downlinks, the SpaceX video stream from a rocket), that's not negligible, when you move at high speeds.
And really, 200 km/h is pretty high speed, for example, for typical cellular comms (phones!). Say, your phone is using a band at 1.8 GHz, and you're sitting in a car going 100 km/h = 28 m/s. Then, your Doppler shift is 1.8·10⁹ Hz · 28/3 · 10⁻⁸ = 168 Hz. That doesn't sound so much, but it means that it takes but 1.5 ms until your phase has rotated by 90°, if my in-head calculations aren't off.
Since high-rate communications usually needs the phase to be consistent, you need to measure your channel way more often than that – so, more than 1000 times a second, your channel changes due to Doppler, if you don't correct the Doppler shift first!
And that, at a relatively benign 100 km/h (LTE does 300 km/h of bullet trains), at a relatively benign 1.8 GHz (modern WiFi works at 5.8 GHz).

How about communications with high velocity spacecraft and probes on other planets?

Of course. When you want to receive what was transmitted by a low-earth orbit satellite (the lower a satellite is, the faster it has to orbit around earth), you first have to correct for Doppler, before you can do much else.


Strange results with all pass compensation

UPDATE: As Hilmar pointed out, OP did convolution in the time domain (as is clear in this code). This answer is misleading so will be deleted).
You would not convolve the two responses as you are in the frequency domain. It is the time domain impulse response of the equalization filter that you would convolve with, and convolution in the time domain is the product in the frequency domain. For the product of complex signals, the phases would add, so in your simple example (which may not be realizable due to causality) the phases would add to zero.
For the product of complex signals, the phases would add:
$$K_1e^{j\phi_1}K_2e^{j\phi_2} = K_1K_2e^{j(\phi_1 + \phi_2)}$$


Overlapping power-of-cosine windows with constant weight

your result is interesting but you should take a look at :

F. J. Harris, "On the use of windows for harmonic analysis with the discrete Fourier transform," in Proceedings of the IEEE, vol. 66, no. 1, pp. 51-83, Jan. 1978.

doi: 10.1109/PROC.1978.10837

Abstract: This paper makes available a concise review of data windows and their affect on the detection of harmonic signals in the presence of broad-band noise, and in the presence of nearby strong harmonic interference. We also call attention to a number of common errors in the application of windows when used with the fast Fourier transform. This paper includes a comprehensive catalog of data windows along with their significant performance parameters from which the different windows can be compared. Finally, an example demonstrates the use and value of windows to resolve closely spaced harmonic signals characterized by large differences in amplitude.
keywords: {Discrete Fourier transforms;Fourier transforms;Frequency;Harmonic analysis;Oceans;Parameter estimation;Signal processing;Signal resolution;Signal sampling;Smoothing methods},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1455106&isnumber=31261

There a number of Figures of Merits (FOMs) that are discussed for windows including $cos^k(),\;
 k=\{1,2,3,4\}$ but not at the generalized overlaps you show.  Looking at the main table in Harris's paper, as $k$ increases, many of the FOMs like the ENBW get worse.  Also while the constant energy is preserved for your result, the coherent gain of the DFT decreases.  The only FOM that touches on overlap is Overlap Correlation.

If you may permit me being blunt, but $cos^{20}$ as windows go, doesn't seem to offer much advantage in general terms. I think that your result is interesting and someone might publish it, but you might hold off going down to the Porsche dealership.


Channel impulse response in freqency domain

You are right. If you have MATLAB, you can do a little experiment to double check yourself:

h = randn(100, 1);     % make some random impulse response
F = dftmtx(length(h)); % make the DFT matrix (100-by-100)
abs(fft(h) - F*h)      % this number should be small


Tied rank adjustment for Spearman Correlation

After finding the two rank vectors Rank_a and Rank_b, just calculate the correlation coefficient of ranks:

c = cov(Rank_a,Rank_b)/(std(Rank_a)*std(Rank_a))

From Wikipedia, you can only use
$$C=1−\frac{6\Sigma(R_x−R_y)^2}{n(n^2−1)}$$ if the ranks are distinct integer values (in this case are not).

  If ties are present in the data set, this equation yields incorrect results.


In addition to what others have already said, I'll try to answer it from a purely practical point of view (this is also a variant of the overlap-add technique).

If your FFT length is 2048, then an overlap of 1024 (50%) means that you have twice as many analysis (FFT) frames (as compared to the number of frames without any overlapping). A 512 samples overlap (75%) means 4 times as
many frames and so on.

The purpose of overlapping is primarily to reduce the effect of windowing. Most windowing functions (e.g Blackman, Hamming etc) are taper-shaped, which means that they drop to 0 (or close to 0) near the frame edges. This, of course, affects FFT results and
we may lose some important information (e.g. transients).

So to reduce this negative effect of the windowing, we use overlapping. The basic idea here is that we can average FFT results form overlapping frames and thus obtain a better frequency representation of
our time-domain signal. The actual frequency resolution is still the same as without overlapping.

As an example, let's say we use 1024 samples for the FFT and we have a 50% overlap (512 samples). We can average 3 overlapping frames to compute the final (averaged) magnitude values

M_avg[frame][i] = (M[frame-1][i]+M[frame][i]+M[frame+1]]i])/3 ; // i is just a bin index--in this case the bins range from 0 to 512

So, to compute the magnitudes for, say, frame #5, we average magnitudes from frames #4, #5# and #6 and divide by 3 (the number of averaged frames)

There are other use case for overlapping but the one described above is probably used the most.


Add multiple aborption coefficients into virtual impulse response

A few points:

This is very poor model to approximate a real room impulse response, so you you should make sure that it can meet your needs before putting work into it. I think this one here may be already better: https://www.mathworks.com/matlabcentral/fileexchange/25965-fast-simulation-of-acoustic-room-impulse-responses--image-source-method-
The notion of a frequency independent reverb time (or absorption) is seriously flawed. In order to get something remotely realistic, you need to bake in some frequency dependency. Brute force would to be calculate one RIR for each octave, band-pass filter and summing  them up.
If you want to expand this code, you should probably look at equation 9.

A quick hack could like like this

c = sqrt(r(1)*r(2)).^abs(e) .*
    sqrt(r(3)*r(4)).^abs(f) .*
    sqrt(r(5)*r(6)).^abs(g)

where $r$ is a vector with the reflection coefficients for each wall matching the x,y,z coefficients. That will be "close enough" for higher order reflection but pretty wrong for lower order reflections, especially if the reflections of opposite walls are fairly different.


Short answer: We can, but it is not always as robust as frequency domain techniques.

Longer answer: This is very broad topic, but let me try to throw some light. Time domain techniques tend to be naive. For example some are based on the RMS energy monitoring of the signal being recorded by an accelerometer. Main disadvantage is that on early stages of fatigue you can't distinguish between faulty and working bearing. Additionally setting a threshold is a real pain. Also some people are using methods based on crest factor (ratio of peak to RMS), knowing that crest factor for impulsive sounds is high. Although there are some 'fancy' methods for analysis in time domain, such as: New Time Domain Method for the Detection of Roller Bearing Defects

Frequency domain techniques (such as Cepstrum or Bispectrum) on the other hand are more robust and have properties that are desirable, such as you can very easily correlate detected fault frequency with a specific bearing and type of damage. For example when you use methods such as Bispectrum, which is (simply speaking) based on detecting presence of higher harmonics. Because it is also taking phase into consideration, in theory even if harmonic is almost covered with noise, you can still detect if its correlated with fundamental frequency - smart isn't it?


Modulation scheme for a noisy channel with a complicated transfer function

Far field, so there is potentially significant reverb

That's basically a multipath channel with a large delay spread.

In communications engineering, the measures taken to revert that effect would be called equalization.

So, that's less of a modulation problem, but a channel estimation problem. Estimating channels with large delay spread / high frequency-selectivity is relatively hard, and so is equalizing them once you know the channel. You can deal with that by sending a known sequence ("preamble") against which the receiver correlates.

Since the result of that correlation is the autocorrelation function of the preamble, convolved with the channel impulse response, you can get pretty far with that. Especially considering you're free to pick a preamble that is spectrally white, i.e. has a sharp (dirac) impulse as ACF.

Based on that, you could build a zero-forcing equalizer, which is essentially simply taking that impulse response, inverting its Fourier transform, and transforming back to time domain.

Problem with zero-forcing equalizers is that you get strong noise amplification at the frequencies the channel is essentially "deaf". In audio engineering, that's what happens when the guy at the FOH doesn't pay attention, and suddenly you get a really unpleasant tone at some frequency when someone starts using their microphone.

Other options include MMSE equalizers, LMS adaptive filters, etc.; there's a large corpus of literature on both equalization for communications, and equalization for making audio work.

I'd actually, for simplicity reasons, recommend a different approach:

Instead of trying to equalize this very long, and hence very frequency-selective channel, you'd simply take it, and divide it into many narrower channels.

If you choose to pick these equidistantly and use a synchronous boxcar-shaped pulse shaping on each of these, you end up with OFDM (Orthogonal Frequency Division Multiplexing).

The idea is that instead of this complex equalizer that you need to use the full channel with a high bandwidth data signal, you get a lot of narrower channels that you divide your data to, and because they are so narrow, an equalizer is just a phase and amplitude correction – much, much easier!

And to boost the benefits, OFDM is mathematically just an IDFT - DFT pair, and can hence be implemented really efficiently using FFTs.

  Typical residential background noise

yeah, so noise. Some frequencies might be especially disturbed, but most of this noise will be relatively broadband.

So, channel coding with a relatively long code might be a good idea. Probably with an outer and an inner code, to convert burst noise to something benign and then throw a high-rate code on the rest.

  phone microphone quality, phone case, and room transfer function can be all over the place

Yeah, so that's covered by equalization.

  microphone handling noise and hand partially blocking the mic

burst Noise and fading.

  There may be mild time-variance, but probably not a lot

well, that's why you'll need to re-estimate your channel regularly.

If you want something to play with: GNU Radio is a software framework for live signal processing and can be directly attached to microphone sources or speaker sinks.

It comes with OFDM blocks, but there's also projects like gr-adapt that offer adaptive filters to compensate channels.

I'd start with the OFDM examples, and just replace the radio hardware block with a block that interpolates the complex baseband to match your audio sampling rate, and shift it up to a carrier frequency before discarding the imaginary part.

Do the reverse (add zero imaginary part, shift down to baseband, decimate to baseband bandwidth) at the receiver  – Tadah, OFDM audio testbed.

Pick an OFDM FFT length that is large enough to divide the channel into subchannels of  a couple dozen Hz bandwidth, at first (that system will then have very long symbols, and thus be sensitive to channel changes, but it should be robust against complex channel impulse responses).

Start with a really "easy" modulation – BPSK. If that works, and a constellation diagram says your SNR permits it, you can go for higher modulations like QPSK.

My guess is that staying within the confines of PSK is probably a good idea. OFDM isn't really easy on the analog chain – it requires a large dynamic range. In fact, I think I remember that if you want to avoid more than one in thousand signals clipping, you need a 10 dB headroom over average power.


You can create your own custom kernel filter using something similar to this example.

if you just want to find image gradients there are other options such as sobel and laplace

If your aim is edge detection, I find canny is best for this in most cases.


Calculate a DFT of a simple finite signal {+1, −1, +1, . . .})

Computing the DFT sum

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

directly can be tedious. A computer would generally use a Fast Fourier Transform (FFT) algorithm.

For the special class of exponential signals

$$x[n]=a^n,\qquad n=0,1,\ldots,N-1\tag{2}$$

$-$ and your example is a special case of Eq. $(2)$ $-$ you can find a closed-form expression for $X[k]$ by using the geometric sum formula. (Just be careful with those values of $k$ for which the denominator of that formula would become zero.)


Derivation of Circular Mean Square Error

The first term in $e^2_p$ is just the variation of the estimates around the true value.

The second term is due to the problem that happens at the end of the periodic region. Sometimes, the noise is enough to move the estimate from $-\pi+\alpha$ to $\pi - \beta$.

For example, the orange plot in the figure below represents the variation of the estimates for a true value around zero.  This variation is captured by the first term.

The blue plot in the figure below represents the variation of the estimates for the true value around the end of the period.  The trouble is that some estimates end up all the way at the other end of the periodic region.  This causes a bias to appear.  The second term in (36) is aimed at capturing that bias.


multiplying the received signal by carrier, in OFDM system

That is just a simulation of a direct conversion downmixer.

You're probably looking at a simulation of OFDM in e.g. Matlab. In a real OFDM system, the signal at carrier frequency is analog, and has to be mixed down to complex baseband with an analog mixer. Of course, in an analog mixer, there's no discrete time $t$, but just an oscillator of the right frequency with which the I and Q output branches of the mixer are generated.

Since complex baseband is so fundamentally important to how OFDM works, I strongly recommend refreshing your knowledge of it and of direct conversion mixers, if this sounds new to you.


function y = rect(x);

    x = abs(x);

    if x == 0.5
        y = 0.5;
    else
        if  x < 0.5
            y = 1;
        else
            y = 0;
        end
    end

end

put this in a file and name it "rect.m".  make sure it's on your MATLAB search path.


Update

I would like to rephrase the original problem with slightly different notation and claim that there is no need to measure the coherence of the product of the circulant matrix A and the DCT matrix B as asked by the OP.

In the compressive sensing setting, the coherence of a sensing matrix $A$ is useful when we are using it for constructing measurements $y = A \alpha$ where $\alpha$ is a sparse vector. Then, the coherence properties of $A$ can be used to give an estimate of number of measurements required in $y$ so that $\alpha$ can be recovered stably from $y$. $\alpha$ is actually a sparse representation of some signal $x$ in some orthonormal basis $\Psi$ which is represented by the equation $x = \Psi \alpha$. In this particular case, $x$ is the image 32x32, $\Psi$ is the IDCT2 transform, $\alpha$ is the sparse representation of $x$ obtained by DCT2 transformation which is the equation $\alpha = \Psi^T x$. Thus, the matrix $B$ from OP is $\Psi^T$. So the relationship between the measurements $y$ and the original signal $x$ is
$y = A \Psi^T x$. The sparse recovery of $x$ from $y$ would happen in two steps.

Construct $\alpha$ from $y$ by using some algorithm $\Delta$. We represent this algorithmic recovery step as:  $\alpha = \Delta(A, y)$. $\Delta$ could be basis pursuit or OMP or anything else.
Construct $x$ from $\alpha$ by using the equation $x = \Psi \alpha$.

In this overall recovery process, the matrix $A B = A \Psi^T$ plays no role. Hence, there is no need to measure the coherence of $AB$.

The old answer was incorrect. The coherence of a matrix $A$ doesn't change if it is pre-multiplied by an orthogonal matrix $B$ as $BA$. Coherence of a matrix can change in arbitrary ways if it is post-multiplied by an orthogonal matrix as $AB$. I am not aware of any particular results of relating the coherence of $A$ with $AB$.

Below is old answer.

You need to compute the coherence of the matrix A. B is an orthonormal matrix. So it doesn't impact coherence of the combined matrix AB .  For computing the coherence of A, normalize its columns and then take pairwise inner products of the columns. The coherence is the highest inner product in magnitude of the pairs of normalized columns of A.

For more information about coherence, see my tutorial notes here .

You may also want to learn about Babel function which is a generalization of coherence.


Fast Fourier Transform showing the same results but with different files?

I'm not familiar with this type of code, but I can take a guess at what you want to do. Your variable names are inappropriate. I suggest you change  your code to something like the following (comments are inside brackets):

spectrum = fft(data)   [Assumes you've transformed a data sequence whose length is 'N' samples.]

fs = input("What's the 'data' sampling rate (in Hz): ")

delta_freq = fs/N   ['delta_freq' is the frequency spacing (measured in Hz)
between the spectral samples contained in the 'spectrum' sequence]

freq_axis = np.linspace(0, delta_freq, (N-1)*delta_freq)

plt.plot(freq_axis[:N//2], np.abs(spectrum)[:N//2]*(1/N))  [Plots the normalized (i.e., the 1/N scaling) spectral magnitude over the positive-frequency range]


Fastest Algorithm for Distance Transform

Pedro F. Felzenszwalb and Daniel P. Huttenlocher have published their implementation for the distance transform. You cannot use it for volumetric images, but maybe you can extend it to support 3d data. I have only used it as a black box.


Here is the poor man's approach:

var histogram = []
foreach pixel
    histogram.append(abs(pixel.value - pixel.south.value))
    histogram.append(abs(pixel.value - pixel.east.value))
sort(histogram)
sharpness = histogram[histogram.count * 0.95]
upscaled = sharpness < X // You choose X experimentally


How does google automatically create 3D cityscapes?

They basically use 3D reconstruction based on ground and sometimes aerial footage of places, generally using stereo cameras. With this footage, they build a 3D point cloud based on the images to do the reconstruction (which is particularly trivial for a stereo setup).

I was told that if the reconstruction isn't good enough for more important/well known landmarks, they will have someone create a better 3D model in some CAD software. But I do think most of what they have, impressively, is done algorithmically based on 3D point clouds and potentially some other things I am unaware of (homographies maybe?).


A receiver needs to synchronize on those two timings in order to demodulate with best performance:

Transmitter symbol rate
Transmitter carrier

Synchronization state of the art way is to use a PLL. You need two loops processing at the same time. The same PLL architecture can be used for both synchronizations but instantiated separately. That is the case for the filtering part of the loop.
Except:

Each loop has a different timing error detector: One detects symbol rate error and the other carrier phase error. Different algorithms for each.
The feedback correction is applied in different areas in the chain: Symbol rate feedback is applied on a resampler, carrier error correction where the frequency translation is done.

Also important, a PLL is a fine synchronization. Which means that if you have a too large offset it won't be able to synchronize alone. A process ahead in the chain should perform the coarse synchronization.
Timing error detector for each loop should be smartly selected. As they can interfere between themselves. For example Gardner when applicable is able to estimate the symbol rate timing error with a carrier phase offset inside.


You mix up the convolution theorem: You cannot multiply the impulse responses, but merely their spectra, if you want to go this way.

However, a simpler way (esp. when you have Matlab available and look for a numeric solution) is to calculate the overall impulse response. For example, considering the filter $x[n] \rightarrow x_1[n]$ you have:

Let $x[n]=\delta[n]$. Then, $x_1[n]=h_1[2n]$, i.e. the impulse response of the overall filter is $h_1[2n]$. Note, that n is zero-based here (i.e. first sample is at $n=0$).

You can also do this in Matlab. Generate a delta impulse and send it through the filter bank. The resulting impulse responses are your effective FIR coefficients for each filter path:

N = 50;
h0 = ones(N, 1).*exp(-(1:N)');            % a lowpass
h1 = sin(3*(1:N)).*exp(-0.3*(1:N));     % a highpass

x = [1; zeros(1,1)];  % a delta impulse

x1 = downsample(conv(x, h1), 2);

x2p = downsample(conv(x, h0), 2);
x2 = downsample(conv(x2p, h1), 2);

x3p = downsample(conv(x2p, h0), 2);
x3 = downsample(conv(x3p, h1), 2);

x4p = downsample(conv(x3p, h0), 2);
x4 = x4p;

% x4 contains the impulse response of the filter x->x4
f = 1/pi*linspace(-pi, pi, 1000);
f1 = 1/pi*linspace(-pi/2, pi/2, 1000);
f2 = 1/pi*linspace(-pi/4, pi/4, 1000);
f3 = 1/pi*linspace(-pi/8, pi/8, 1000);
f4 = 1/pi*linspace(-pi/8, pi/8, 1000);

subplot(2,1,1);
plot(f, mag2db(fftshift(abs(fft([h0, h1'], 1000, 1)), 1)));
xlim([0, 1]);
xlabel('f/pi');

subplot(2,1,2);
hold off;
plot(f1, mag2db(fftshift(abs(fft(x1, 1000)))));
hold on;
plot(f2, mag2db(fftshift(abs(fft(x2, 1000)))), 'r');
plot(f3, mag2db(fftshift(abs(fft(x3, 1000)))), 'k');
plot(f4, mag2db(fftshift(abs(fft(x4, 1000)))), 'g');
hold off;
grid on;
xlabel('f/pi');
xlim([0, 1]);

As you see that the sample rate is reduced for each path, the maximum frequency is also reduced by half. If you were to plot all frequency responses over the full frequency range, you'd need to periodically repeat them.


OpenCV Haar Classifier - best type of positive images

You probably want to train on as tightly cropped rectangular region around the bee as possible.

I'm not sure if you already thought of this but you should remember that this classifier can only handle recognizing objects in rectangles whose axes are aligned with image, i.e. no rotated rectangles. This is going to be tricky if the bees can be facing any which way. A human face detector has the advantage that faces are almost always found right side up in the image. If you wanted to recognize upside down faces too you would probably train only on right side up images and then check your input image twice--once with the original input and once with it flipped upside down.

You probably want to do something similar for the bees. All of your training images should be aligned so that the bee has approximately same orientation and scale. Then if you want to detect different orientations you can rotate your input image.


Matlab filter design,through group delay

The Solution finally was to match the input vectors in the group delay simulating Matlab functions.That means that the frequencies and the group delay Vectors given as input in Matlab should be SHORT and not as long as i had to give,and exactly on the same length (which makes sense).The most tricky part is to find which length was that would fit in and run the simulation.


Are radar resolution and detection capabilities not very tightly related?

Resolution is (usually) referred to as the ability to distinguish two closely separated returns, not targets.

The distinction is how far along the processing chain of the radar you are, resolution, given by the inverse of the bandwidth is often limited by the hardware (waveform generation, filters, amplifiers, sampling..) and the environment. While Radar detection is done as a further processing to distinguish between what you want to see (targets) and clutter (birds, ground, sea, trees, ..) edit: and is usually limited by the signal to noise/clutter ratio. The detection step will usually contain multiple looks.

You can therefore have say a really good resolution, but since all you see is "noise/clutter", you can't really detect the target you are looking for. Edit: Conversely you may have a really good detection algorithm that uses say excessive averaging and multiple 'looks', but if your resolution is poor, two closely spaced targets may appear as a single target.

Edit: Note that I use noise and clutter interchangeably, as it simply refers to what you don't want to see.

Edit: To simplify things: resolution will depend on the bandwidth, while detection will depend on noise level and processing.

I have attempted to sketch a generic radar system below, the resolution is given by the receiver output, while the detection algorithm will typically use multiple looks (for averaging and doppler discrimination) and use a threshold to alert the user of one or more target.


Audio EQ Cookbook without frequency warping

here is a quick look at how the 5 degrees of freedom for the parametric EQ can be viewed.  it's my take on what Knud Christensen of tc electronic came up with about a decade ago at an AES convention and this patent.

so, forget about the Cookbook (and the issues of Q and bandwidth therein) and consider (in the s-plane) the parametric EQ as the sum of a bandpass filter (with a $Q$ value) in parallel with a wire:

$$ H(s) = (G_\text{boost} - 1)\frac{\frac{1}{Q}\frac{s}{\omega_0}}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} + 1 $$

$G_\text{boost} = 10^{\frac{dB}{20}}$ is the gain of the peak (or valley, if $dB<0$).  the gain at DC and at Nyquist is 0 dB. that's a 2nd-order IIR and there are 3 independent parameters. 2 more to go.  so we next add an overall gain parameter:

$$ H(s) = G_0 \left((G_\text{boost} - 1)\frac{\frac{1}{Q}\frac{s}{\omega_0}}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} + 1 \right) $$

that's 4 knobs to twist.  one more knob to add (without raising the filter order) and we will be done with adding more independent parameters.

so what Knud does here is replace that "wire" (that trailing "$1$" in the transfer function) with a prototype shelving filter that must have the same poles, the same $Q$ and $\omega_0$ as the BPF, so that the denominator is the same.  the transfer function of that shelf is:

$$ H_\text{shelf}(s) = \frac{R\left(\frac{s}{\omega_0}\right)^2 + \frac{\sqrt{R}}{Q}\frac{s}{\omega_0} + 1}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} $$

where $R \triangleq 10^{\frac{tilt}{20}}$ and $tilt$ is the gain differential of the shelf in dB.  this is what offsets the gain at Nyquist to be different than the gain at DC.  After bilinear transformation, Nyquist gets boosted by $tilt$ dB and the gain at DC remains unchanged.  Like the $dB$ boost parameter, the $tilt$ parameter can be either positive or negative.  $G_0 \cdot R$ is the linear gain at Nyquist.

put that all together and you get:

$$ \begin{align}
H(s) & = G_0 \left((G_\text{boost} - 1)\frac{\frac{1}{Q}\frac{s}{\omega_0}}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} + H_\text{shelf}(s) \right) \\
 & = G_0 \left((G_\text{boost} - 1)\frac{\frac{1}{Q}\frac{s}{\omega_0}}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} + \frac{R\left(\frac{s}{\omega_0}\right)^2 + \frac{\sqrt{R}}{Q}\frac{s}{\omega_0} + 1}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} \right) \\
& = G_0 \frac{R\left(\frac{s}{\omega_0}\right)^2 + (G_\text{boost}+\sqrt{R}-1)\frac{1}{Q}\frac{s}{\omega_0} + 1}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1}  \\
\end{align}$$

no matter how you look at that, this has 5 degrees of freedom and those 5 biquad coefficients are fully defined from these 5 parameters.  doesn't matter if you map from $s$ to $z$ using the blinear transform or the trapezoidal rule (effectively the same thing) or any other method that does not change the order of the filter.  you may have to fudge the definition of $Q$ or bandwidth, you may have to compensate $\omega_0$ and/or $Q$ for frequency warping effects (like you get with the blinear transform), but if you paid big bucks for something that gets you a 2nd-order IIR filter, it doesn't matter if you implement it with any Direct Form or transposed Direct Form or Lattice or Normalized Ladder or Hal Chamberlin's State-Variable or Andrew Simpson's modeling of linear analog with trapezoidal integration, eventually you get to 5 coefficients and they can be mapped to these 5 independent parameters.  it's all the same.  whether or not you paid money for a license or not.  the math is stronger than any claims made by whomever you are licensing from.

Just FYI, i solved where the true peak or valley frequency is when there is a $tilt$ that is not zero.  the frequency where the peak or valley has been nudged over by the tilt is:

$$ \omega_\text{peak} \ = \ \omega_0 \
\sqrt{ \frac{Q^2 \left(R - \frac{1}{R}\right)}{G_\text{boost}^2 - R + 2 Q^2 (R - 1)} + \\ \sqrt{\frac{G_\text{boost}^2 - \frac{1}{R} + 2 Q^2 \left(\frac{1}{R} - 1\right)}{G_\text{boost}^2 - R + 2 Q^2 (R - 1)} + \frac{Q^4 \left(R - \frac{1}{R}\right)^2}{\left(G_\text{boost}^2 - R + 2 Q^2 (R - 1)\right)^2} } } $$

you can see that when $tilt = 0$, then $R=1$ and consequently $\omega_\text{peak} = \omega_0$.   the peak gain $G_\text{boost}$ might also have to be adjusted a little and that has yet to be worked out.  a good first guess would be $G_\text{boost} \leftarrow \frac{G_\text{boost}}{\sqrt{R}}$ or maybe $G_\text{boost} \leftarrow G_\text{boost}-(\sqrt{R}-1)$.


best meathod to lock onto a cdma signal?

CDMA signals typically have a pilot that is sent along with the data.  This pilot is a known chip sequence that is transmitter over and over again.  Receivers can use the pilot to lock onto the signal by correlating the received signal with the expected pilot.

The autocorrelation is typically done by FFT'ing the received signal, multiplying it by the FFT'ed expected pilot, and IFFT'ing the result.  This is done both for computational efficiency and because the circular correlation that you get with FFT's is actually helpful in this case because the pilot is repeated.  Thus you should be able to "catch" the pilot no matter which phase it starts at in your captured signal.

Once you have locked onto the pilot you know the signal's phase and it should be pretty straightforward at that point to despread the signal and recover the data, providing you know the data's spreading sequence.


Deriving Time-Series Data from PSD when the Inverse Fourier Transform doesn't converge

Even if you had succeeded in deriving the inverse Fourier transform of $1/|f|$ (note that you need the absolute value) then you still wouldn't have a time-series with a $1/f$ behavior. What you would get is the auto-correlation sequence.

A good way to generate pink ($1/f$) noise is to generate (approximately) white noise with a random number generator, and filter that noise with a filter whose frequency response approximates $1/\sqrt{f}$ over a large frequency range. Of course the frequency response needs to have a finite value at $f=0$, but it is possible to design a filter approximating $1/\sqrt{f}$ for $f>f_{min}$ with a relatively small value for $f_{min}$. If you have a question as to how to design such a filter you should probably formulate this as a new question. You can also have a look at this site.

And by the way, the inverse Fourier transform of $1/|f|$ does exist in terms of distributions: check this.


Your understanding of the convolution process is correct. However, note that in the convolution the kernel is first mirrored before the dot product. Though, for a symmetric kernel this does not matter.

Understanding the convolution in your way, of finding specific features, is a similar interpretation as that of a Matched filter. Furthermore, using the convolution to find a specific part of the image (e.g. eye), is sometimes called Template matching (but here, usually correlation is used, which is equal to convolution except the kernel is not mirrored).

Another interpretation of the convolution is done via the convolution theorem, which says that in the frequency domain, the frequency response of the filter and the image are multiplied. Essentially, you can then understand the filter to enhance / reduce higher or lower frequencies, (or do more fancy stuff in the frequency domain). With this interpretation, it sometimes becomes more understandable why a certain kernel has a certain effect (e.g. the Gaussian filter is easier to understand as a lowpass than a filter to look for Gaussian shapes in the image).


Viewpoint dependant deformation of objects in photos

Since you have listed "illumination", I am assuming that "condition changes" include changes in lighting.  So I would add cast shadows, which will change if you move the light source.

Also, all the geometric transformations you have listed are affine. But if you move the camera perspective changes as well. So you may want to consider the projective transformation.


Using I/Q demodulation/LIA for a fixed-frequency application of several signals - do the phases of the input signals affect the result?

Under the assumption that your words "extract the amplitude of those four components" actually mean "extract the spectral magnitudes of those four components", then the initial phases of the original four periodic signals do not matter. It's straightforward to prove this behavior by modeling your down-conversion process using software.


Unknown symbol/expression in text about adaptive filters (cst)

The standard normalized step-size LMS algorithm computes the current step-size according to

$$ \mu = \frac{c}{s_k^T \cdot s_k}  $$

where $c$ is a suitable scale factor and $s_k^T \cdot s_k$ is the total energy of the current tap inputs. The algorithm aims to adjust step size according to input signal power; when input has large power then decrease the step-size so that filter coefficient update will not deivate much. But if the signal goes to zero or very small levels (silence durations, or small variance sections), then the division will possibly be very large and algorithm may produce unbounded (very large) outputs or even diverge...

To prevent this from happening, a small nonzero factor; the regularisation factor, is added to the denominator

$$ \mu = \frac{c}{s_k^T \cdot s_k + \delta} $$

to limit the maximum of the fraction to $c/\delta$ when input power goes to zero.


That sure looks like your are including the .WAV file header as part of your data.

When you read your file in, and write it out, using the file operations, you are not accounting for the header.  The header can vary in size (it's actually three RIFF headers), but is usually 44 bytes long.

These are the headers in my own (C++) words:

struct FirstHeader
{
char     RiffID[4];
int      RiffLength;
char     WaveID[4];
char     FormatID[4];
int      FormatLength;
};

struct SecondHeader
{
short    Always0x01;
short    TrackCount;
int      SamplesPerSecond;
int      BytesPerSecond;
short    BytesPerSample;
short    BitsPerReading;
char     Filler[16];
};

struct ThirdHeader
{
char     DataID[4];
int      DataLength;
};

Do a search on "Wav file header" and you can get lots of details.  If you don't need any of the parameters, I suppose just looking for the "data" in the header and the next 4 bytes are the UINT32 value of the length.  So your data starts the byte after.

Otherwise, follow the RiffLength and FormatLength to get to the DataID field.

The file can be mono or stereo, 8 16 or 24 bit, etc.


How to design an inverse for a Butterworth low-pass filter?

The simplest and probably best approach is to use a second order IIR filter exchanging numerator and denominator from the original filter - as you suggested - but move the double pole at $z=-1$ (Nyquist) slightly towards $z=0$, i.e. slightly decrease the pole radius to make the filter stable. Just play around with the exact value. The higher the value, the closer to Nyquist will the approximation of the inverse filter be.

This little Matlab/Octave program shows how to do it:

[b,a] = butter(2,.005);     % some 2nd order Butterworth LP filter
r = 0.98;                   % desired pole radius < 1
b2 = a;                     % EQ numerator
a2 = [1,2*r,r*r];           % EQ denominator
k = sum(a2)/sum(b2);
b2 = k*b2;                  % scaling
[H1,w] = freqz(b,a,512);    % Butterworth frequency response
H2 = freqz(b2,a2,512);      % EQ frequency response

% plot EQ'd response
plot(w/pi,20*log10(abs(H1.*H2)))
grid on, xlabel('normalized frequency'),
ylabel('dB'), title('equalized response')


What can be done with I/Q data which would otherwise be impossible?

The main motivation for using quadrature (I/Q) signals is spectral efficiency.
Most signals are naturally baseband: they have no carrier, their energy is concentrated around DC, and they are limited to bandwidth $B$. An example would be an audio signal, which has bandwidth of roughly 20 kilohertz.
To transmit such a signal, often it must be modulated (or "upconverted") using a carrier. For example, in AM broadcasting, the audio signal is first limited to a bandwidth $B=10\,\text{kHz}$, and then upconverted to a frequency of roughly one megahertz.
The issue is that, when upconverted, a signal occupies a bandwidth $2B$, which is double the original bandwidth $B$. So, an AM signal occupies $20\,\text{kHz}$ "on the air". The extra 10 kilohertz do not have any additional information; they are simply a loss of half of the available bandwidth.
Since bandwidth is a very limited and expensive resource, such a waste is unacceptable.
One way around this issue is to use single side-band modulation (SSB), which eliminates the waste. However, SSB signals are more difficult to generate and have other disadvantages.
The most popular solution is to use quadrature signals. This allows transmitting two signals of bandwidth $B$ over a chanel of width $2B$, restoring the average spectral occupancy to $B$ hertz per signal.
To see more of the math behind quadrature transmitters and receivers, see this answer.


Simulating room impulse response between rooms

This is an FDTD based, wave equation solution. It is the best free solution I have found. None of the other answers has helped so far.


Find the valley in the given in a sequence of numbers

The following code will find the lowest point between two prominent peaks.

%Smoothing the data
movAvgFiltOrder = 8;
b = filter(ones(1,movAvgFiltOrder)/movAvgFiltOrder,1,sums);

% finding peaks from the smoothed curve
[vals,pos] = findpeaks(b);
idx = find(vals>mean([min(vals),max(vals)])); % removing peaks which are below the average,
%this is helpfull if some lower amplitude peaks are caught at sides of the prominent peaks.

peak_pos1 = pos(idx(1)); % farthest prominent peak on one side of valley
peak_pos2 = pos(idx(end)); % farthest prominent peak on other side of valley

% correcting the peak positions for filter delay (Coarse tuning)
peak_pos1 = peak_pos1 -  movAvgFiltOrder/2;
peak_pos2 = peak_pos2 -  movAvgFiltOrder/2;

% Looking for any other peaks nearby (Fine Tuning)
tuning_range = min(movAvgFiltOrder,floor(abs((peak_pos2-peak_pos1-1)/2)));
p1 = peak_pos1-tuning_range:peak_pos1+tuning_range;
p2 = peak_pos2-tuning_range:peak_pos2+tuning_range;
[~,r1] = max(sums(p1));
[~,r2] = max(sums(p2));
peak_pos1 = p1(r1);
peak_pos2 = p2(r2);

% Identifying valley
valley = sums(peak_pos1:peak_pos2);

% finding minimum point in valley
[min_val,pos_min_val] = min(valley);
valley_amplitude = min_val; % amplitude of minimum point
valley_position  = peak_pos1+pos_min_val-1; % finding the valley point in the entire array

This valley_position will give the poisition of the lowest point between major peaks.


Given a system with Transfer Function and its desired output. Is it possible to find the required Input?

The transfer function is given as ...

Your transfer functions seems wrong. It's not "proper", i.e. there ought to be a $1$ in the denominator. You can try to rewrite it this way, but it would make it non-causal. Not sure if that's a typo or an error of your System Identification

  However, the inversed is not proper

No, your original transfer function is improper, the inverse is just fine.

  If it is possible, could you kindly give me some broad instruction

Make sure your original transfer function is proper, fix this first
Check if your corrected transfer function is causal, stable & minimum phase, i.e. all poles and zeros inside the unit circle. If that's true, you can simply invert it by exchanging poles and zeros.
If you have zeros outside of the unit circle, you can still invert this, but the result will be non-causal. You can still approximate an inverse by either adding bulk delay (if your application can tolerate that) or by least square error time domain solution.


Modeling channel response as complex vs splitting I and Q

After giving it more thought, I see my original confusion. We need to be careful about what we mean by the channel's impulse response. Normally for an equalizer for this application, we are not modeling the channel alone, but rather also including the frequency mixing up and down, before and after, transmission through the channel. When we include the mixers in this way, we get the cross coupling as I mentioned for the first case. So if $h(t)$ is the response of the channel alone, then $g(t)=e^{-j\omega_ct}h(t)$ is the response of the entire system from baseband to baseband. The second case with independent responses would be for applications where the channels are not mixed together at any stage and are transmitted and received independently.


What is shown in the plot is the PSD (Power Spectral density), which is measured in W/Hz (or dBc). So, in order to get the actual power in a given frequency range, you'd need to multiply with the bandwidth of each bin (which is determined by the length of the DFT window). Let the bandwidth of each bin be $B$.

Then, your formula becomes

$$
THD = 10\log\left(\frac{P_{HD}}{P_S}\right)=10\log\left(\frac{\sum_i B\cdot S_i}{B S_{HD}}\right)
$$

where $S_i$ denotes the Density at the $i$th harmonic and $S_HD$ denotes the density at the fundamental. As you see, B cancels in the expression, so you can ignore it. Note, that $S$ are in linear scale, (i.e. you'd need to transform the measurements in the diagram to linear scale before summing them up etc.)


maximum output of a LTI system within linear limits

For a given LTI system, …, after which it becomes non linear

never. If you're modeling a system to be LTI, it's LTI: Linear (and time-invariant).

You can of course model a system to be non-linear in general, but linear for a range  of input. But then the answer to your question is "that's in your own model, and we can't tell you"!


Can PSNR > 48 used to validate the signals between two audio files?

If the reconstructed audio signal is exactly same as original signal then MSE =0. And if Max pixel value is 255 (8-bit representation), then the value of PSNR = 20*log(255) = 48dB

That statement is plain wrong. If the two signals are identical, then MSE = 0, then division by zero is undefinied and you don't get a PSNR at all, but "infinity".

  Since the PSNR of the songs are more than 48, can I conclude that the decoded audio signal in the wav file is valid?

PSNR is really not a very useful metric; you can get these values with a codec that simply broken, or one that sounds great.

Let's throw a bit of math in here: let's say we map our 8 bit to - to +1. Our original audio signal $S$ is well-scaled and its amplitude is pretty normally distributed¹, centered around 0, and has a variance $\sigma^2=\frac14$, so that clipping very rarely happens, i.e. $S\sim\mathcal N(0;\frac14)$.

Then, what's the PSNR of a codec that simply always produces "0"?

\begin{align}
\text{PSNR} &= 20\log_{10}\left(\frac{\text{Maxval}}{\sqrt{\text{MSE}}}\right)\\
&=20\log_{10}\left(\frac{\text{1}}{\sqrt{\text{MSE}}}\right)\\
&=10\log_{10}\left(\frac{\text{1}}{\sqrt{\text{MSE}}}\right)^2\\
&=10\log_{10}\left(\frac{\text{1}}{\text{MSE}}\right)\\
\end{align}

So, what's our MSE? It is

\begin{align}
\mathbb E\left((0-S)^2\right)&= \mathbb E\left( S^2\right)\\
&S \text{ is zero-mean}\\
&=\text{Var}\{S\}=\sigma^2\\
&=\frac14
\end{align}

Therefore,

\begin{align}
\text{PSNR} &= 10\log_{10}(4) \\
&\approx 6
\end{align}

So, what this would prove is that your codec is better than just producing the mean value of the information-theoretically worst possible continuous audio source – but that's about it.

¹ that assumption about the amplitude distribution is slightly wrong, but wrong in favor of your codec.


Unit Problem in Designing a Filter for a Given Auto Correlation Function

Radians are considered to be dimensionless.  See Are angles dimensionless? and Dimensionless quantity.  They are considered to be pure numbers like pi.

So $\alpha$ is in Hz, which is a measure of 1/second, and $s$ is also considered to be measured per second.


Modelling an existing DSP filter by analysing the differences between its inputs and outputs

If I understood correctly, you have a signal $x[n]$ and an unknown discrete-time LTI filter, so that you can look at the output of the filter $y[n]$. Now you are looking for the impulse response of the filter $h[n]$.

The output follows a convolution rule
$$y[n]=x[n]*h[n]=\sum_{m=-\infty}^{\infty}x[m]h[n-m]$$
The process to find $h[n]$, given $x[n]$ and $y[n]$ is called deconvolution. There are some approaches to do that.

For example, if the Fourier transform of the signals exist, using convolution theorem, the convolution can be converted to product in the frequency domain: $$Y(\omega)=X(\omega)H(\omega)$$
Therefore, $$H(\omega)=\frac{Y(\omega)}{X(\omega)}$$
You should however be careful about the zeros of $X(\omega)$. If $x[n]$ can be an option, then you should choose it such that it can excite the system at all frequencies such that you can probe the desired frequency range that you are interested in $H(\omega)$. Alternatively, you can excite the filter in repeated number of times at different frequencies, which can be interpreted as sampling the frequency response. This will give you a sampled spectrum, so if you know the filter does not have singularities between the sample frequencies, then you can model it by an LTI filter with your desired order (depending on your acceptable approximation error).

In response to your second question, after you acquired an approximation such as $\hat{h}[n]$ for the impulse response, then you can apply it to any desired input by using convolution.


transmitting the real/imaginary parts of pass-band signal

Will we be able to recover the complex modulated symbol based on $s_2(t)$ too?

Yes, no information is lost between the two.

How can we prove that mathematically?

One way to think about it is to write out each individually and see what you have. Let $a(t)=a_I(t)+ja_Q(t)$, expanding $s_1(t)$ and $s_2(t)$ we have:
\begin{align}
s_1(t) &= a_I(t)\text{cos}(2\pi f_ct)-a_Q(t)\text{sin}(2\pi f_ct) \\
s_2(t) &= a_I(t)\text{sin}(2\pi f_ct)+a_Q(t)\text{cos}(2\pi f_ct)
\end{align}
The in-phase and quadrature parts are present in each so we just need to bring them back to baseband to get $a(t)$ back. Further, whatever workings you have to do this for $s_1(t)$ you can reuse for the processing of $s_2(t)$ by first applying a simple $+\frac{\pi}{2}$ to it ($\text{cos}(\theta + \frac{\pi}{2})=-\text{sin}(\theta)$).


Matlab FFT keeps plotting mirrored values

It is a mathematical property of the DFT that its output is conjugate symmetric if its input is real-valued. That means that for the magnitude of the DFT output signal $X_k$ of length $N$ it holds
$$
|X_{N-k}| = |X_k|.
$$
So what you see is to be expected. (In your example $3200\,\mathrm{Hz}/2 - 266\,\mathrm{Hz}\approx1330\,\mathrm{Hz})$ If you're only interested in the positive frequencies, just plot the first $N/2$ elements of the DFT output, the rest is redundant.


Finding the maximum and minimum frequency of a 2D image

I think our friends on Stackoverflow have a solution for you that would be helpful. Despite this being MATLAB code, I believe the solution extends to a general answer to your question.


Some time ago I wrote some codes just to test some filters methods, and i could see how can be hard make certain frequencies disappear completely, some methods cause ripple effects in the signal and when you apply FFT on the filtred signal to see the frequency response, you note that the frequency you want to eliminate still be there, attenuated, but you can still see!

Lets go, three(3)  sinusoidal concatenated at 250 hz, 400 hz and 900hz sampled in 8000hz and I will try to eliminate 900hz from the original signal!

My First test, AKA a very simple equalizer in Frequency Domain, calculate the dB attenuation 10.^(dB/20) and multiply by the FFT bin that you want eliminate, for dB attenuation -24dB my result was:

My Second test EQ biquad filter by Robert Bristow-Johnson (time domain), for dB attenuation -24dB my result was:

My Third test biquad LPF by Robert Bristow-Johnson (time domain), my result was:

And my last test using Sinc Filter (Frequency Domain), my result was:

And that's all, in short, everything depends on what you really need !


How to solve equations for the MMSE receiver?

I doubt there is general closed-form solution. For the first equation, you essentially want to solve $$\mu=\sum_k \frac{1}{a_k+b_k \mu}$$ for $\mu$ which equals to search for the zeros of a $k+1$th order polynomial, which is impossible for $k>3$. So, I fear you need to do a numeric solution. Once you have $\mu$ you can go for a simple calculation of $\kappa$.

Regarding a numerical solution, here's a plot of $$\mu-\sum_k \frac{1}{a_k+b_k \mu}$$ which should become $0$ for the solution:

N = 7
d = np.random.rand(N)
K = N+1
M = 16
p = 1

f = lambda mu: mu - 1/(K-1)*sum(1/(M*p*d_i*(1-(K-1)/M+(K-1)/M*mu)+1) for d_i in d)

mu = np.linspace(-5,5,10000)
plt.grid(True)

plt.plot(mu, f(mu))
plt.ylim((-0.5,0.5))

So, there are multiple solutions. Maybe you need to pick to positive one?


How do I convert a real baseband signal to a complex baseband signal?

To convert a real signal sampled at rate $2B$ to its complex baseband representation (sampled at rate $B$), you want to map the frequency content in the range $[0, B)$ in the real signal to the range $[-\frac{B}{2}, \frac{B}{2})$ in the resulting complex signal. This can be done in a couple different ways:

Design a linear filter to approximate a Hilbert transform. Run your real signal $r[n]$ through the filter to yield the transformed signal $\tilde{r}[n]$. Use this to form the analytic signal:
$$
r_a[n] = r[n] +j\tilde{r}[n]
$$
$r_a[n]$ will contain only the positive frequency components of the original real signal $r[n]$; all of the negative frequencies will be zero (assuming a perfect Hilbert transformer; in practice, the effect will not be perfect). Thus, you have isolated the desired frequency band $[0, B)$.

Multiply $r_a[n]$ by $e^{-j\frac{\pi}{2}n}$ to effect a frequency shift of $-\frac{B}{2}$, shifting the desired frequency content to the range $[-\frac{B}{2}, \frac{B}{2})$. Then, decimate the signal by 2 by discarding every other sample. The result is a complex baseband signal sampled at rate $B$.
My preferred approach is a more straightforward implementation of the shift that you're looking for:

Multiply $r[n]$ by $e^{-j\frac{\pi}{2}n}$ to effect a frequency shift of $-\frac{B}{2}$, shifting the desired frequency content to the range $[-\frac{B}{2}, \frac{B}{2})$. The result is a complex signal that is centered in the appropriate place (the center of the band is at zero frequency).
Apply a lowpass filter to pass only the content in the range $[-\frac{B}{2}, \frac{B}{2})$ (to the extent required while meeting your application's antialiasing requirements).
Decimate the signal by 2 by discarding every other sample. The result is a complex signal sampled at rate $B$.

In practice, I always use some variant of strategy #2. You can make it even more computationally efficient by implementing the decimation as part of the lowpass filtering process, for instance with a polyphase filter.


What's the physical meanning and principle of applying anisotropic diffusion in transform domain？

Physically, there is usually not much implied meaning (but this may depend on a case-by-case basis). Transforms like shearlets are actually used to embed the data into anisotropic or edge-aware transform domains.

Mathematically, it may make more sense. Isotropic diffusion tends to create blur (smear) in whatever domain it is applied. Thus, if you want to obtain a result that is compact in the considered domain, be it transformed or not, applying non-isotropic diffusion is very tempting.

Another way to see it mathematically, is to understand that isotropic diffusion corresponds to solving a least square problem that will create many small non-zero coefficients by trying to spread the error everywhere. However, a helpful and widely available assumption about signal is their sparsity in the transform domain, i.e., they have few non-zero coefficients in a suitably chosen transform domain that concentrate the energy. This corresponds to solving problems where the exponent of the norm of the prior term is less than 2 (typically 1) and thus the use of anisotropic diffusion to find a solution.


How to separate mixed signals, 3 identical gears meshing together

Yes it would make sense to analyse the sound of the vibrations without the gear interactions if you can, and then delete that sound profile from the overall sound using a noise remover that works with a sample of the noise only prior to the noise + recorded signal. You should say roughly what frequencies you expect from the gears based on the RPM and number of teeth, because it affects analysis methods. On this and many other programs there is "noise reduction using noise sample profiling" http://wiki.audacityteam.org/wiki/Noise_Reduction it is tremendously efficient.

There are no obvious correlations or patterns in the graph, it is true. Read about engine vibration analysis. I found a paper on engine sound graphs when i was researchign SPWVD, and indeed it is 100ds of times more precise than FFT and pseudo smoothed wigner ville distribution makes waterfall graphs which show the clearest patterns of sound that are physically measurable.

To seriously analyse your vibrations, you should not use coarse FFT, instead use pseudo smoothed wigner ville distribution which is available as an addon. this program can do it too, see the graphs: http://www.christoph-lauer.de/sonogram


What is the relationship between scales and frequency in a Morlet wavelet transform?

I found some details, for those who might come after me.

Scale and frequency should be inversely proportional, which is as one might figure. However, there seems to be an issue (possibly a bug) in PyWavelets that makes this not always the case.

There is an internal function in the PyWavelets library called scale2frequency() that takes 2 arguments plus an optional 3rd: wavelet, scale, and then precision. By its name, I assume that it's made to take a scale and return the associated frequency. However it's an involutory function: $f(f(x)) = x$. So it can also be used to calculate scale from frequency. (This might be undocumented behavior, and I have no idea if it'll still work this way in future versions.)

The problem I was running into involved the optional precision argument. It would appear that as the wavelet's central frequency increases, the required precision became very high. I was able to mitigate the problem somewhat by increasing precision beyond the default value of 8. But as it gets higher and higher, the amount of memory and processing power increases drastically, to the point where I burned through all my system's virtual memory long before it was high enough to return valid answers.

My final workaround involved altering my data's sampling rate by a factor of 1000 or so. In essence, I just told it that my data was in the Hz range rather than the kHz range. This made it so that the code wasn't returning any obvious errors. It's not ideal, but it seems to work well enough.


When should I use a frequency spectrum rather than a spectrogram representation?

The spectrogram suffers from a time/resolution restriction - the smaller your time step, the smaller your analysis windows become, the worse your frequency resolution gets.

It is a waste of frequency resolution to break down a stationary signal into tiny temporal slices. Thus, if you know in advance that your signal you analyze is stationary; or if you have pre-segmented your signal into stationary segments (using a transient/onset detector; novelty detection; or an auxiliary signal/sensor), there is no reason to lose frequency resolution by using the spectrogram.


How to undo the quantization step in DCT?

From the code it looks like the quantization step is doing two things: 1) dividing by "Q", and 2) rounding.  Undoing the divide by Q is easy.  The rounding is where the actual quantizing is happening though, and unfortunately that is the part that you cannot undo- at least not without more information.  Rounding causes information to be lost and you cannot undo that.