What are the advantages of Laplace Transform vs Fourier Transform in signal theory?

Laplace transforms can capture the transient behaviors of systems.  Fourier transforms only capture the steady state behavior.  Of course, Laplace transforms also require you to think in complex frequency spaces, which can be a bit awkward, and operate using algebraic formula rather than simply numbers.

If you want to see the power distribution of a signal over time, Fourier transformations are often the easiest way to do it.  However, if you want to understand what a system does when you flip a light switch, you typically need Laplace transforms.


In order to get the actual frequency of each bin, not normalized frequency, you must compute the following:

fres = fs/fft_size; % Bin spacing
f = fres * (0:fft_size/2); % Bin centers and frequency axis


Difference between delta modulation and differential PCM

Since the encoded values are differences between original values, the values are relatively small values and they are in general distributed as small positive and negative values around zero. This helps to further compress the data by using some algorithm that can use shorter codes to represent smaller values that appear frequently in the data, and longer codes for larger values that appear less often.


This can be accomplished with a convolution of the distance transform.

Use a distance transform on the edge of the mask. Then threshold this distance transform to remove values beyond some distance. I think the secret to getting the shading is to convolve the distance transform result with a kernel that looks something like this:

[ -1.0  -1.0  -1.0
  -1.0   0.0   0.0
  -1.0   0.0   1.0 ]

This should get you started in the right direction:

#include "opencv/cv.h"
#include "opencv/highgui.h"

using namespace cv;
using namespace std;

int main() {
    Mat mask, dist, bevel;
    mask = Mat::zeros(200, 400, CV_8U);
    rectangle(mask, Point(30,30), Point(180,180), Scalar(255), -1);
    circle(mask, Point(30,30), 50, Scalar(0), -1);
    circle(mask, Point(180,180), 50, Scalar(0), -1);
    circle(mask, Point(300,100), 75, Scalar(255), -1);
    imshow("1",mask);

    //find edges and invert image for distance transform
    Canny(mask, dist, 50, 150);
    dist = 255-dist;
    distanceTransform(dist, dist, CV_DIST_L2, CV_DIST_MASK_5);
    threshold(dist, dist, 20, 20, CV_THRESH_TRUNC);
    blur(dist, dist, Size(3,3));
    dist.convertTo(bevel, CV_8U);
    equalizeHist(bevel, bevel);
    imshow("2",bevel);

    //convolve with secret sauce
    float d[] = {-1,-2,-3,
                 -2, 0, 0,
                 -3, 0, 1 };
    Mat kernel(3, 3, CV_32F, d);
    kernel = kernel - mean(kernel)[0];
    filter2D(dist, dist, CV_32F, kernel);

    //normalize filtering result to [-1, 1]
    double maxVal;
    minMaxLoc(dist, NULL, &maxVal);
    dist = 128 * dist / maxVal;

    //convert and display result
    dist.convertTo(bevel, CV_8U, 1, 128);
    bevel = bevel.mul(mask)/255;
    imshow("3", bevel);

    waitKey(0);
}


Squared magnitude of System Function H(z)

It's important to understand that the equation
$$|H(z)|^2=H(z)H\left(\frac{1}{z}\right)\tag{1}$$
is only valid on the unit circle, i.e., for $z=e^{j\omega}$. For complex-valued systems the general form of $(1)$ is
$$|H(z)|^2=H(z)H^*\left(\frac{1}{z^*}\right)\tag{2}$$
For real-valued systems we have $H(z)=H^*(z^*)$ and, consequently, $(1)$ is valid for real-valued systems. Note that both $(1)$ and $(2)$ are only valid on the unit circle.
The reason why we use $(1)$ or $(2)$ as a generalization of $|H(e^{j\omega})|^2$ to the whole complex plane is that we want an analytic function that is actually a $\mathcal{Z}$-transform of some sequence. And the function $H(z)H^*(z)$ does not satisfy that requirement.
Note that the IDFT of $H^*(e^{j\omega})$ is $h^*[-n]$. So if we find the $\mathcal{Z}$-transform of $h^*[-n]$ we have what we need:
$$\sum_nh^*[-n]z^{-n}=\left(\sum_nh[n](z^*)^{n}\right)^*=H^*\left(\frac{1}{z^*}\right)\tag{3}$$


I'll assume that you want to do online monitoring.

You can always store the coming data from the GUI in a vector in MATLABA.
Then you can run simple code to plot them.

This also depends on the type of measurements you have. If the set is capable on detecting power ranges as well you should have the following readings:

High Alpha, Low Alpha, High Beta, Low Beta, High Gamma, Low Gamma, Delta, Theta, Attention, and Meditation.


How to triangulate 2 sources on a 2d plane formed by 4 sensors?

Conceptual partial answer:

You can use 3 mics to find position in 3D space, except you can't tell if it's above the plane or below the plane (unless you already know there's nothing above the plane, of course).

Since your 4 microphones are in the same plane, the 4th microphone is redundant and doesn't help you know whether the source is above or below the plane.  (Though it improves SNR/accuracy of estimate, as would any number of additional mics in the same plane).

There is a delay due to speed of sound from the source ($S$) to each microphone.  The delay difference between each pair of microphones gives you a hyperboloid (green cross-section) that the source must be on, in order for the waves to be delayed by that relative amount (due to extra distance $a$).  Each pair of mics defines a hyperboloid, so by combining 3 pairs of mics, you can find the intersection of the 3 hyperboloids, which will define a unique set of 2 points.  (Symmetrical above and below the plane.)

(4 mics gives you 6 pairs of microphones, so by averaging the intersection point of all 6 hyperboloids you can improve the accuracy of the estimate, but you still don't know if it's above or below the plane.   With a tetrahedron of mics you would only have one intersection point.  You could probably weight the average, too, since some source positions will be more accurately measured by some pairs of mics, etc.)

Also, the sound pressure drops off as 1/r with distance.  The attenuation (amplitude) ratio between each pair of microphones gives you a sphere (yellow cross-section) that the source must be on.  By combining this information from every pair of microphones you can find average intersection of all 6 spheres and improve accuracy from that, too.

The spheres and hyperboloids are usually not tangent at the intersection point, so they each improve accuracy in a particular direction:

To find the relative delay, you can use cross-correlation of the signal from the 2 mics.  There will be a peak at the "lag" at which the two signals line up with each other.

To find the projection of the source location on the plane itself, just throw away the z coordinate.

See https://en.wikipedia.org/wiki/Multilateration#Solution_algorithms for more details on actual implementation.


multiplication of a function with a Fourier-transformed equals to Fourier-transformed with a function

Hint:
You only need to use the definition of the Fourier transform and its inverse transform to show that if
$$G(f)=\mathscr{F}\{g(t)\}$$
then
$$g(f)=\mathscr{F}\{G(-t)\}$$
holds. Then just use the result from $b$.


MRI K-space to image: why track frequencies in two dimensions?

in MRI, you don't start with the image as pixels in spatial domain, but with the k-space data!
So, you don't "construct k-space from the image", you "construct the image from k-space", which takes a 2D inverse Fourier transform.
And that's not "do it row-by-row"; that is "do it on all rows, then all columns of the result".


Why is cepstrum used in MFCC instead of autocorrelation sequence?

Usually in speech signal processing cepstrum is used to represent low and high frequency components , which are multiplied with each other(in time domain its a slowly varying signal convolved with high frequency signal).

  What's the reasoning behind using log value?

logarithm is taken to remove the multiplicative effect and convert it to addition of two components. $$\log(a\cdot b) = \log(a) + \log(b)$$

  Why does that make a difference?

Now the two components are in addition form.

  Why isn't the autocorrelation sequence used in MFCC, LPCC, PLPCC, etc.. instead of the cepstrum?

Same reason, because we need to remove multiplicative effect of slowly varying component and fast varying component.

In speech signal processing slow varying component represents different phonemes and thus act as feature vectors for many applications(ex ASR). hence we need to remove the effect of fast varying component to isolate slow varying component.


Wavelet transform: How to compute the initial coefficients when only samples are available?

I remember that using the discrete samples was termed  "wavelet crime" by Gilbert Strang (Strang and Nguyen, Wavelets and Filter banks, 1996, pages 232 sq.), or "wavelet sin" by others.

On prefiltering discrete signals to obtain "better" wavelet coefficients, here are some references I have in store:

Wavelet coefficient computation with optimal prefiltering, IEEE Trans. Signal Processing, 1994, Xia, X.-G. and Kuo, C.-C.J. and Zhang, Z.
On the initialization of the discrete wavelet transform algorithm, IEEE Signal Processing Letters, 1994, Abry and Flandrin
From the wavelet series to the discrete wavelet transform - the initialization, IEEE Trans. Signal Processing, 1996, Xiao-Ping Zhang and Li-Sheng Tian and Ying-Ning Peng
Sard-optimal prefilters for the fast wavelet transform, Numer. Algorithms, 1997, S. Ehrich
A new prefilter design for discrete multiwavelet transforms, IEEE Trans. Signal Processing, 1998, Xia, X.-G.
Quadrature prefilters for the discrete wavelet transform, IEEE Trans. Signal Processing, 2000, Johnson, B. R. and Kinsey, J. L.
Meaningful MRA initialization for discrete time series, Signal Processing, 2000, D. Veitch and M. S. Taqqu and P. Abry
Efficient Wavelet Prefilters with Optimal Time-Shifts, IEEE Trans. Signal Processing, 2003, Ericsson, S. and Grip, N.

It seems however that no prefiltering does not induce too much errors for compression, denoising, deconvolution. But for detection, or using multiple  different  wavelets at once (like with dual-tree wavelets, or unions of wavelet bases), or with multiwavelets, it is quite mandatory.


Nyquist Theorem adding two same frequency near to Nyquist Frequency with phase shift

When you add 2 or more sinusoids at the same frequency $f_o$ but with different phase shifts you get a sinusoid at same frequency but an additional attenuation term. Mathematically, you will have following:
$$cos(2\pi f_ot + \phi)+cos(2\pi f_ot + \theta) = 2cos(\frac{\phi - \theta}{2}).cos(\frac{2\pi f_ot + \phi + 2\pi f_ot + \theta}{2})$$$$= 2cos(\frac{\phi - \theta}{2}).cos(2\pi f_ot + \frac{\phi + \theta}{2}),$$You see basically, you have an attenuation term $2cos(\frac{\phi-\theta}{2})$ which will change the amplitude of the sum of the sinusoids, but the frequency of the sinusoid remains at $f_o$. Hence, if you were sampling it with $f_s>2f_o$, then even after adding them you will not have any aliasing.

Only thing of concern is that because of attenuation term driving the amplitude of the sinusoidal, your signal might become 0 now, as in you example of $2^{nd}$ sinusoidal at the same frequency $f_o$ but having phase difference($\phi - \theta$) of $\pi$ compared to $1^{st}$ signal. That means $2cos(\frac{\pi}{2})$ which is equal to 0, and hence you have actually nothing to sample, both signals are completely out of phase and cancelling each other.


Weighted Nuclear Norm Minimization for Image Denoising

Most of the Denoisers in Image Processing make a simple assumption - The data has small number of freedom degrees while noise has high number.
Hence if we try to represent the given noisy data with small number of parameters we probably match the data while exclude most of the noise.

In the case above we try to limit the number of Singular Values of the base which represents the patches of the image.
The assumption means that most of the patches can be represented by a dictionary with small number of Singular Values -> Small number of degrees of freedom to create "Crazy Patches".


The DTFT is continuous, so you have to approximate the DTFT on a computer.

Zero padding “interpolates” the DFT, so  a lot of zero padding will approximate the DTFT


Do FIR filter have any application in control theory?

FIR filters are fairly common in some areas of control theory. As they usually incur a lot of added phase/time-delay, they are not really usable in the feedback path of regular control systems, but they are useful when the added phase/time-delay is not affecting the system in an adverse way, or when the particular phase response and time-delay is desired.

Examples:

Feed-forward control. FIR filters are useful for producing filters that approximate arbitrary frequency responses, hence they can be used to shape a reference signal. A typical example is to use an FIR filter with the inverse frequency response of the plant -- trying to counteract the dynamics of the plant in order to get a desired output. Phase/time-delay is not interfering with the stability or performance since the computation can be done offline. FIR filters can often produce higher performance than IIR filters, especially where there are non-minimum phase zeros.
Learning systems (iterative learning control [ILC] and repetitive control [RC]). In these systems, again, FIR filters are often used to generate an inverse response filter of the plant. This is done to improve the performance of the learning, as this enables a higher gain to be used for a larger bandwidth. The time-delay is here desired (RC), or it does not matter, as the computations are done offline (ILC). Low-pass FIR filters are also used here for their linear-phase properties, as the linear-phase property makes the filter not interfere with the phase of other filters, and this improves performance.


What are integration scale and differentiation scale?

The $\sigma_I$ determines the scale level at which the Harris corners are computed. Coarser scales (higher values of $\sigma_I$) correspond to larger corners.

The $\sigma_D$ is effectively the window size, over which the derivatives are summed to generate the entries of the matrix. If $\sigma_D$ is too small, then the detector will be seriously affected by noise.  If it is too large, then small corners may wash out.  In the original non-scale-adapted Harris detector you only have $\sigma_D$.

By the way, this is the first time I see the terms "integration scale" and "differentiation scale", and I am not sure I like them. The "differentiation scale" is actually the Gaussian envelope, over which the derivatives are integrated. :)


why is moving towards higher frequencies so important in frequency spectrum?

So, if frequency increases, signals possesses higher energy and can
  travel far.

Not really, it's more the other way round:

When transmitting signals with radiowaves in free-space, most telecommunications engineers use the Friis equation which gives the amount of received power $(P_r)$ according to emitted power $(P_t)$, distance between emitter and receiver ($R$) and wavelength ($\lambda$) :
$$
\frac{P_r}{P_t} = G_tG_r\left(\frac{\lambda}{4\pi R}\right)^2
$$
$G_t$ and $G_r$ being the transmitter and receiver's antenna gain.
We see with this equation that $\frac{P_r}{P_t}$ is proportionnal to $\frac{1}{f^2}$ (because $\lambda = \frac{c}{f}$) which means that the higher the frequency, the higher the loss (and ^2 higher !).
The equation $E = h\nu$ is more related to quantum physics and has no practical use in telecommunications. Maybe it means that higher frequency signals are "more energetic" that low frequency ones, I don't really know honestly, maybe someone else* can answer better.

So the need to increase frequencies is mostly due to :

Saturation of the currently used spectrum : if you read international and national regulations and standards, you'll see that being able to emit on the RF spectrum is very expensive, because it's actually saturated, every frequency is reserved for a certain application (e.g. GSM signals is my country are transmitted in 200 KHz wide canals around 900 Mhz-ish frequencies). Still in my country, a new mobile phone operator had to pay billions of € to the governement to access some mobile comm. bands. Increasing the frequencies will give everyone more "room" to transmit radiowaves.
Increasing data rates, because having a higher carrier frequency also allows to widen the communication channel
Miniaturization (higher frequencies = lower wavelengths = smaller antennas and systems)

*edit :

  your claim on signal energy and frequency relation, which is based on photo-electricity (E=hv) effect is improper. First of all, the physical nature of waves has nothing to do in developing the mathematical notion of waves and their spectral representation by periodic signals. And even though it's true that the magnitude-squarred |x(t)|2 signal power relation is basicly adopted from power calculations in physics, that power (and energy) is implied to be the energy carried out by the waveform of a classical physics kind.

from Fat32's comment above


How to plot the Hilbert Spectrum in Hilbert-Huang transform?

Final step is pretty straightforward. All you need to do is to apply the Hilbert Transform to each IMF and extract the instantaneous frequency from analytical signal. Instantaneous frequency is given by:
$$\omega(t)=\dfrac{d\phi(t)}{dt} $$
where $\phi(t)=\mathrm{arg}[x_a(t)]$ (unwrapped phase of the analytical signal). Keep in mind that MATLAB (Octave) hilbert function already returns analytical signal.

In fact you might want to re-use this code from Mathworks File Exchange: plot_hht it is a good starting point.


Is it possible to discretely sample the function

yes, with spacing between samples closer than $\frac{\pi}{\alpha}$.
no.  it's not sufficiently "bandlimited"
yes, the convolution of a sinusoid with anything is a sinusoid of the same frequency.  so it's "yes" for the same reason as 1.


Sharpen Defocused Image (Deconvolution / Image Restoration)

You're after an algorithm in the family of "DeConvolution".
Specifically in your case, is called Blind Deconvolution.

Yet if you have some assumption the Blur you can use Wiener Filter or Lucy Richardson.
Both of them are actually the MMSE Estimator just with different assumption of the noise.
Both of the methods are actually "Inverse Filter" on an Low Pass Filter, which means they are High Pass Filter, just like the Unsharp Filter you applied.
The difference is those methods are optimal per given "Blur Model".

Yet, if you go after the Blind Deconvolution, today the best estimators are based on some prior of the image.
Something like the distribution of the Gradient of the image and stuff like that.

Yet if you are only after the text in the image.
Something simple like Wiener Filter + Edge Detection should do the work to give you most data.


How to determine the central frequency from a vectorial spectrum analyzer?

The S11 represents the relative amount of signal that is reflected from the antenna, so the lower the better, because lower means more radiated by the antenna.

Usually, we consider the antenna to be "working" at a given frequency if its S11 is lower than -10 dB (of course, this requirement can be changed according to the application. spatial antennas are more around -20 dB for example).

Therefore, your antenna seems to be working at three distinct frequencies (with the s11 < -10 dB criterion). Is this a problem ? It depends on the hardware you put behind it during operation.


A sequence OFDM modulated with normalization factor

The power of the non-zero elements is $|1+j|^2=2$, so the average power is

$$\overline{P}=\frac{12\cdot 2}{52}=\frac{24}{52}=\frac{6}{13}$$

So in order to get unity average power, you have to normalize by $\sqrt{13/6}$.


what is nyquist rate of $h(t)\cdot h(t)$ and $h(t)\circledast h(t)$

squaring a signal effectively doubles the bandwidth of the signal (and doubles the Nyquist rate, which is the lower limit of the sampling rate).  convolving a signal with anything, even itself, does not in general increase the bandwidth.


Refer good books for the following topics in Image processing

for Multiresolution based fusion discrete wavelet transform you may refer the following
" A Wavelet based image fusion tutorial " by Gonzalo Pajares ( available at: www.ELSEVIERcomputer science.com)

"Wavelet based image fusion techniques - An introduction ,review and comparison by Krista Amolins

for Curvelet transform refer " Wavelet for image fusion" by shih-Gu Huang


$x[k]$ is the input signal to the filter; it could be anything. There is no "rule" for what signal you can input to a linear filter. With that said, the stability of filters is usually talked about in a bounded-in, bounded-out (BIBO) sense. If the input signal is unbounded, then it's likely that the output signal will be also. So if anything, the "rule" that you're looking for would be to only apply signals $x[k]$ whose magnitude $|x[k]|$ is bounded below some value $B < \infty$.


Treating ADC quantization noise as a uniform white noise source is generally valid.  Please refer to this post where I further detail practical limitations on using white noise assumptions for quantization noise as well as oversampling considerations: What are advantages of having higher sampling rate of a signal?
It is typically good design practice to place the ADC quantization noise to be at least 10 dB below the amplified thermal noise of the input signal to minimize the ADC quantization impact on overall noise figure (As the waveform we wish to measure/capture typically is the input waveform including its signal and noise components, not the ADC quantization noise!).
When choosing 10 dB for example as I noted above, the ADC contribution to increasing the noise floor (and hence the noise figure given the same signal level) using the white noise assumption for quantization noise is 0.4 dB given as:
$$10 \log_{10}(1+10^{-\text{dB}/10})$$
where dB is 10 in this case. (This is simply from adding the relative powers and converting back to dB).
The trade off here is dynamic range vs noise figure degradation. Certainly you can place the quantization noise to be much lower below the input amplified noise floor to minimize overall increase in noise (by providing more gain prior to the ADC) but that would be at the expense of dynamic range. For example, if instead of placing the quantization noise 10 dB below the input amplified noise floor as we did above, we placed it 30 dB below, we would lose 20 dB of dynamic range in the process.  Typically I use numbers in the range of 10 dB to 15 dB as starting parameters for input noise floor vs ADC quantization noise, but ultimately is part of an overall noise figure budget.
This is identical to a cascaded noise figure computation for those familiar with that in a receiver design. In this case the ADC has an equivalent noise figure, thus constraining the overall gain requirements prior to the ADC in order to meet an overall noise figure objective while balancing instantaneous dynamic range constraints.
Below are my slides where I walk through a typical ADC receiver design showing how this was used in terms of a "Noise Figure" for the ADC. If there is interest in the rest of this (which shows the same considerations for the maximum input signal), please ask as an additional question specific to that and I will post that there as well (since I already have the slides):

Note: The "Noise Figure" shown in the graphic below refers to the "front-end" noise figure which is the entire receiver up to the input of the ADC. This predicts the amplified noise floor level since that would be the receiver gain up to this point plus the receiver noise figure. This is not the "ADC Noise Figure" that I reference in the text and prior graphics.


Except Logarithmic operation, other operations could be done on combined signal or separate signals both will give same result.

$STFT$($n$) $+$ $STFT$($S$) $=$ $STFT$($n$ $+$ $S$)
$mel$$(STFT$($n$) $+$ $STFT$($S$)) = $mel$($STFT$($n$ $+$ $S$)))

so just removing effect of logarithm and adding should do.
though I am not sure that artifacts will go or not, but you need not do reverse Mel(I am assuming that mel filter banks central frequencies and bandwidths are same for both $S$ and $n$).


Low pass equivalent of bandpass white noise

Note that for the complex noise envelope $z(t)=x(t)+jy(t)$, the autocorrelation $R_z(\tau)$ is defined by (cf. Eq. $(4.1.47)$ in Proakis)
$$R_z(\tau)=\frac12E\big\{z^*(t)z(t+\tau))\big\}=\frac12\big[R_x(\tau)+R_y(\tau)\big]+j\frac12\big[R_{xy}(\tau)-R_{yx}(\tau)\big]\tag{1}$$
As shown in the chapter you refer to, for the real-valued bandpass noise $n(t)$ to be stationary, the following must be true:
$$\begin{align}R_x(\tau)&=R_y(\tau)\tag{2}\\R_{xy}(\tau)&=-R_{yx}(\tau)\tag{3}\end{align}$$
Plugging $(2)$ and $(3)$ into $(1)$ we obtain
$$R_z(\tau)=R_x(\tau)+jR_{xy}(\tau)\tag{4}$$
Note that this expression is different from the one you suggested in your question.
If the power spectrum of $z(t)$ is even, as is the case for bandpass white noise, the autocorrelation function $R_z(\tau)$ must be real-valued, and, consequently, $R_{xy}(\tau)=0$. But note that we don't require $R_y(\tau)=0$.
As a final note, defining a complex envelope is only useful for bandpass processes, such as bandpass white noise, but not for ordinary white noise, because the latter is no bandpass process.


Computationally Efficient Ways of Determining If Pixels Are Clumped

I've not tried this, but I'm wondering if something like a simple low pass filtering, followed by a column sum will do.

No time to explain code right now, but this seems to do something like what you want. Will return in a few hours and explain what it's doing.

R Code Below

library(imager)

N <- 128

noise <- array(runif(N*N*1*1),c(N,N,1,1)) #5x5 pixels, 1 frames, 1 colours. All noise
small_clumps <- as.cimg(noise)
blurry <- isoblur(small_clumps,5)

layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(colSums(blurry-mean(blurry)))
plot(small_clumps)
plot(blurry)

large_clumps <- small_clumps
large_clumps[65:75, 65:75] <- 1
large_clumps[15:25, 35:45] <- 1
blurry_large <- isoblur(large_clumps,5)
# par(mfrow=c(3,1))
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(colSums(blurry_large-mean(blurry_large)))
plot(large_clumps)
plot(blurry_large)


I would try techniques related to low-rank approximations of your data $X$. For instance, the paper GoDec: Randomized Low-rank & Sparse Matrix Decomposition in Noisy Case attempts to decompose
$$X = L + S + G$$
with a low-rank part ($\mathrm{rank}(L) \le r$), a sparse part ($\mathrm{card}(S) \le k$) and a residual $G$.

It has been employed for background modeling, with backgrounds related to the low-rank part, plus sparse en independent objects.

Some video demos are:

GoDec and Robust PCA for background modeling on Escalator
GoDec and Robust PCA for background modeling on Hall

There are alternatives in the field of robust matrix factorization, and I understand it as an instance of sparse coding.


Can bandwidth of the closed loop be bigger than bandwidth of individual elements in it?

Yes, that's perfectly possible and to be expected.

Below is a toy example that starts with a Butterworth low pass filter and puts it in the forward path of a feedback loop. The feedback path has a gain of two.

This makes the whole loop unstable, but it also makes the bandwidth very high.

The original filter's frequency response is below.

And this is the closed loop response.

R Code Below
    #30035

library("signal")

bf <- butter(5, 0.1)
freqz(bf)
dev.copy(png, 'Q30035/Q30035-Butterworth.png')
dev.off()

k <- 2
bf2 <- bf
bf2$a <- bf$a + k*bf$b
bf2$b <- bf$a*k
freqz(bf2)
dev.copy(png, 'Q30035/Q30035-Feedback.png')
dev.off()


How to convert wave from real to complex and vice versa?

Rather than saying the questions you are asking are non-sensical, I will say they reveal a lack of understanding of complex numbers.  There is lots of material to be had with a few simple searches.  I would also recommend that you read my blog article The Exponential Nature of the Complex Unit Circle which gives an explanation of Euler's equation which is absolutely essential to understand in order to comprehend the meaning of DFT bin values.

For your first question:  Technically, a real number is a complex number with a zero imaginary component.  So setting the real part equal to your values and the imaginary part to zero is the correct method.  When you multiplied it by $ i^{0.4} $ it should have rotated each value a tenth of a cycle.  Thus the amplitude stays the same, but the phase shifts by $\pi/5$ radians.

For your second question:  A complex number consists of two values.  A real number consists of one.  Therefore converting a complex number to a real number is a non-sensical proposition.  However, you can calculate the magnitude, which is a real number.  You can also ask what the real part is, which is also a real number.  The imaginary part is a real number multiplied by $i$.

A real tone is a sinusoidal, flat as a pancake.  A complex tone is a spiral, like a slinky.  The definitions I use in my blog articles are:

$$ S_n = A \cos( \alpha n + \phi ) $$

For a real tone, and:

$$ S_n = A e^{i (\alpha n + \phi)} $$

For a complex one.

Notice that a real valued tone is just the average of two complex tones:

$$ \frac{ A e^{i (\alpha n + \phi)} + A e^{-i (\alpha n + \phi)} }{2} = A \frac{ e^{i (\alpha n + \phi)} + e^{-i (\alpha n + \phi)} }{2} = A \cos( \alpha n + \phi )$$

Because

$$ \cos( \theta ) = \frac{ e^{i\theta} + e^{-i\theta} }{2} $$

Hope this helps.

Ced


The blue curve does not seem accurate and is a little non-smooth. However, I will give a single-pole filter to approximate it.

It is not difficult to show that the step response of the filter $$H(z)=\frac{1-p}{1-pz^{-1}}$$
is in the form of

$$s[n]=(1-p^{n+1})u[n]$$

To find the step response, it suffices to shift the signal given in the question to the left and divide it by $2$. Hence, the output at zero is $0.25956786/2= 0.1298$. To find the parameter $p$, substitute it in the given equation for $s[n]$ at $n=0$:

$$0.1298=1-p\Rightarrow p=0.8702$$
So the filter we were looking for is $$H(z)=\frac{0.1298}{1-0.8702z^{-1}}$$
More clearly, $b=0.1298$ and  $a=[1,-0.8702]$.

Here is a comparison in Matlab:

It can be seen that the first $10$ outputs are very close to the actual output. However, after the output becomes non-smooth, there is a small error.
Alternatively, we could use the term at $n=22$ for our reference. In such case, $$1.95277071/2=1-p^{23}\Rightarrow p=0.8497.$$ This will lead to an accurate approximation at large $n$, but an error at small $n$.

Using the code

y = [0.25956786, 0.48537242, 0.68182373, 0.85275161, 1.00148237, 1.1309067, 1.24353671, 1.34155607, 1.42686391, 1.50111103, 1.62198067, 1.67093897, 1.71355379, 1.81104326, 1.83550966, 1.85680735, 1.87534702, 1.89148581, 1.90553486, 1.9177649, 1.92841125, 1.93767929, 1.95277071]
fun = @(P)2*(1-P.^(1:23))-y;
p = lsqnonlin(fun,.5);

$p=0.8644$ gives the least squares (LS) solution to this problem which is compared to the original output in the following plot ($\text{SSE}=0.02796$):

Note that, if the forward filter is $(b,a)$, then $(a,b)$ is the inverse filter which gives you the original input. If the error in the step response cannot be fixed, you can apply the LS approach on the overall output (taking into account the inverse filter) to get a better result.


How to reduce PAPR with constellations choice?

Peak to Average Power Ratio (PAPR) is defined as:

$$
\text{PAPR} = \frac{\text{Peak Power}}{\text{Average Power}}
$$

If you're working with a $M$-PSK constellation, then the peak power is $P_s$ and the average power is also $P_s$ assuming each symbol is equally likely. For $M$-PSK, the $\text{PAPR}=1$.

If you working a $M$-QAM constellation, then the peak power is $\text{max}_{i \in \mathbb{C}} |s_i|^2$, where $\mathbb{C}$ is the constellation. Again assuming all the symbols are equally likely, the average power is $\frac{1}{M}\sum_{i \in \mathbb{C}} |s_i|^2$. The specifics should be worked out for a constellation but is hard to write down because in general there are many possible QAM constellations.

One thing that is directly clear is that $\text{PAPR}_{\text{QAM}}>\text{PAPR}_{\text{PSK}}$. In a QAM constellation, you have some symbols with less power, some symbols with more power, and there is an average somewhere in between. The point is that there are some symbols with more power, and this makes the peak power part of PAPR greater than the average power part. This is in contrast to a PSK constellation, where the symbols lie on a circle and no symbol has more or less power than each other.


Can a 1-dimensional signal be in a domain other than the time or frequency domain?

Really, there's nothing special about time and frequency domain.

The math doesn't care whether you're transforming amplitude over time, gravel over mountain height, or smell intensity over fridge temperature. To the math, you map elements from a function space to elements from a function space when you do the Fourier transform. Full stop.

Any scalar physical entity can be used as basis for a 1D signal. Spatial signals are the most obvious: what's the height profile of something over distance, or the temperature over distance, or color over distance, whatever over distance. In some contexts, the resulting entities are called "spatial frequencies", but sometimes something else. Names are just names.

Or what is the temperature of something over pressure? That's another 1D signal. And it makes a lot of sense to Fourier transform it, if you're dealing with compressible media in a container, things will oscillate, and not necessarily over time. Do I know whether that has a name? No. Did some PhD give it a name in his dissertation about pressure in containers? Pretty likely, if you ask me. Does that name matter to anything but the notation in that dissertation? Maybe. Maybe not.

Very commonly observed: impulse / state space its Fourier "dual" location space in solid state physics.

We can go on with this all day, because, as said anything can be a 1D signal, and some have Fourier transforms with special names, other don't. It's by pure coincidence that you've heard of the time/frequency pair first. If you were a math student, you might have heard about probability densities and characteristic functions first – yet another 1D signal/FT pair.

So, really, to cite someone I think Tukey,

  One can Fourier transform anything – often meaningfully.

Your question is just asking for an arbitrary list of names. Names don't actually matter – they just make it easier to communicate about something. Whether I call something a "time signal" or a "frequency signal", however, makes little difference. For all we care, I could be calling it "original space" and "image space" (under the Fourier transform), because that's what "time" and "frequency" are. It's just convenient that these two have a name.


Question about path loss exponent estimation from measurment

The path loss exponent does not depend on the reference distance. There are two mistakes in your path loss formula. First, path loss does not include the gains of transmit or receive antennas. Second, the path loss exponent must also be part of the expression for the path loss at the reference distance.

According to the Friis transmission equation  the ratio between the received power $P_r$ and the transmitted power $P_t$ is given by

$$\frac{P_r}{P_t}=G_tG_r\left(\frac{\lambda}{4\pi d}\right)^2\tag{1}$$

where $G_t$ and $G_r$ are the gains of the transmitter and receiver antennas, respectively, $\lambda$ is the wavelength, and $d$ is the distance between the transmit and receive antennas. The squared term in parentheses (or its inverse) is the free-space path loss. Note that it does not include the antenna gains. In order to account for other environments than free-space, a path loss exponent $n$ is introduced:

$$\frac{P_r}{P_t}=G_tG_r\left(\frac{\lambda}{4\pi}\right)^2\frac{1}{d^n}\tag{2}$$

In free space we have $n=2$; in lossy environments, $n$ can be greater than $2$.

If we introduce a reference distance $d_0$, Eq. $(2)$ becomes

$$\frac{P_r}{P_t}=G_tG_r\left(\frac{\lambda}{4\pi}\right)^2\frac{1}{d_0^n}\left(\frac{d_0}{d}\right)^n\tag{3}$$

In dB we get

$$\begin{align}10\log_{10}\left(\frac{P_r}{P_t}\right)&=10\log_{10}(G_t)+10\log_{10}(G_r)+\ldots\\&\ldots +20\log_{10}\left(\frac{\lambda}{4\pi}\right)-10n\log_{10}(d_0)-10n\log_{10}\left(\frac{d}{d_0}\right)\end{align}\tag{4}$$

where

$$PL_0=20\log_{10}\left(\frac{\lambda}{4\pi}\right)-10n\log_{10}(d_0)\tag{5}$$

is the path loss at the reference distance $d_0$. Note that in $(4)$ the reference distance $d_0$ will cancel out, so the estimate of the path loss exponent $n$ will not depend on $d_0$. On the other hand, the estimate of $n$ will have a (weak) dependence on the antenna gains. This is however not in conflict with the fact that the path loss as a function of distance is a straight line on a log-log plot.

Note that in practice, the path loss at reference distance $d_0$ is often just set to the free-space path loss, i.e using $n=2$ in Eq. $(5)$ $[1]$. In that case the reference distance $d_0$ will no longer cancel out in $(4)$ and the estimate of the path loss exponent $n$ will depend on the choice of $d_0$. Alternatively, $PL_0$ can also be measured at distance $d_0$, or its can be optimized together with $n$ to obtain an optimal fit to measured data $[1]$.

[1] Wireless Communications, A. Goldsmith, 2005 Cambridge University Press, p. 41.


Your equations are correct (as a first order approximation).

The phon is an attempt to establish a perceptual loudness measurement. Perceived loudness is rather complicated since it needs to take into account what happens in the human auditory system and in the human brain.

Here is is how it works: First you measure the sound energy at the listener's position. The relevant physical quantity is the variation in air pressure, called "sound pressure" and measured in Pascal or Pa like any other pressure. This is then referenced to 20 micro Pascal and turned into dB like this
$$ x = 20\cdot log_{10}\frac{p}{20\mu Pa} $$
where p is the measured sound pressure and x is the sound pressure level in dBSPL. This is a direct measurement of the physical sound energy.

The human auditory system has different sensitivity at different frequencies. It's most sensitive around 2kHz-4kHz, a lot less sensitive at 100 Hz and at 30 kHz you don't hear anything at all, regardless of how much sound energy is there. The "phon" takes this into account. It's based on measurements that compare the perceived loudness of sine waves at different frequencies. The results are often referred to as "equal loudness contours" (see for example http://en.wikipedia.org/wiki/Equal-loudness_contour).

A loudness of X phons means "as loud as a 1 kHz sine wave at X dBSPL". For a 1 kHz sine waves phons and dBSPL are the same thing. For sine waves of other frequencies phons and dBSPL are related through the equal loudness contours. For more complicated signals calculating phons is a bit more tricky.

Phon is not a direct loudness measurement. Perceived loudness, as most human modalities, follows roughly the Weber-Fechner law (e.g. http://en.wikipedia.org/wiki/Weber%E2%80%93Fechner_law). As a first order approximation it takes ten times the energy to double the perceived loudness (all other things being equal). An increase in 10 phone doubles the loudness and an increase in 20 phone does, indeed, increase the perceived loudness by a factor of 4.

To calculate phon (relative or absolute) from a wav file, the file MUST be calibrated in units of Pascal (sound pressure) since the equal loudness contours are non linear and change with absolute level.


Seeking advice on how to denoise a video from poor CCTV lighting

Something pretty basic to do with good generalization would be using one of the following:

The Bilateral Filter.
The Non Local Means Filter.

They are easy to implement and have a decent performance on many models.
Good resource about denoising is - Which Noise Reduction Algorithms Are Used in Commercial RAW Image Processors?


Applying Nyquist theorem to digital vs. analog audio quality

The Nyquist-Shannon reconstruction theorem applies to perfectly bandlimited signals and perfectly flat reconstruction filters.  Unfortunately, neither exists in the real world.  Nor do digital system that use an infinite number of bits per sample to avoid quantization noise.

The headroom allowed by using higher sample rates than 44.1k may allow a closer approximation (for the bandlimit and reconstruction filter response) to the required perfection.  So does using more bits than 16 to reduce quantization noise.

But analog recording is also not perfect, as noiseless analog recording media and head amplifiers do not exist (because of the existence of atoms and non-zero K temperatures).

So there can be a trade-off.  Choose your preferred sounding reproduction artifacts(s).


How does the assumption that symbols are equi-probable hold

If I draw a number uniformly between zero and one, what is the
  probability that they are equal? Mathematically, it should be zero but
  I don't recall why? Can somebody please help in explaining why the
  probability should be zero?

Avoiding formal definitions of (Lebesgue) probability measure, an informal way is thinking the probability at a point of a uniformly distributed random variable is $\frac{1}{\textrm{number of points}}$. Then the number of points of real values between zero and one is infinity, thus the probability must be $\frac{1}{\infty} = 0$.

  For example, there is still a chance,  that when generating six random
  numbers they will be all 1 or all 0? So if the probability is zero
  then why in digital communications do we assume that the symbols are
  $iid$ and that if the symbols are 0 and 1, 50% of them will be 1 and
  the rest will be 0 i.e., why do we assume equi-probable symbol
  probabilities when theoretically they cannot be equal?  In Matlab, I
  tried my best to generate a large sequence of 0/1 but never did I get
  equi-probable symbol probabilities.

The previous arguments change if we divide the interval $[0,1]$ to several sub-interval of width $\Delta$. Thus the probability falling on each sub-interval is $\frac{1}{1/\Delta}=\Delta \neq 0$.

In your example, $\Delta = 0.5$ yields two sub-intervals and each has probability $0.5$. That's where the equiprobable symbol assumption holds.

Experimentally, e.g. with Matlab, you (almost surely) never get exactly the probability $\Delta$ because the number of tests is always finite.


The first step is to verify that both your starting sample rate and your target sample rate are rational numbers.  Since they are integers they are automatically rational numbers.  If one of them wasn't a rational number it would still be possible to make the sample rate change, but it is a much different process and more difficult.

The next step is to factor the two sample rates.  The starting sample rate, in this case, is 44100, which factors to $2^2*3^2*5^2*7^2$.  The target sample rate, 16000, factors to $2^7*5^3$.  Thus, to convert from the starting sample rate to the target rate we must decimate by $3^2*7^2$ and interpolate by $2^5*5$.

The previous steps have to be done no matter how you want to resample the data.  Now let's talk about how to do it with FFT's.  The trick to resampling with FFT's is to pick FFT lengths that make everything work out nicely.  That means picking an FFT length that is a multiple of the decimation rate (441, in this case).  For the sake of the example, let's pick an FFT length of 441, though we could have picked 882, or 1323, or any other positive multiple of 441.

To understand how this works it helps to visualize it.  You start out with an audio signal that looks, in the frequency domain, something like the figure below.

When you're done with your processing you want to decrease the sample rate to 16 kHz, but you want as little distortion as possible.  In other words, you simply want to keep everything from the picture above from -8 kHz to +8 kHz and drop everything else.  That results in the picture below.

Please note that the sample rates are not to scale, they are just there to illustrate the concepts.

The beauty of picking an FFT length that is a multiple of the decimation factor is that you can resample simply by dropping portions of the FFT result, and then inverse FFT what is left.  In the case of our example you FFT 441 samples of data, which gets you 441 complex samples in the frequency domain.  We want to decimate by 441 and interpolate by 160 ($2^5*5$), so we keep the 160 samples that represent the frequencies from -8 kHz to +8 kHz.  We then inverse FFT those samples and presto!  You have 160 time domain samples that are sampled at 16 kHz.

As you might suspect, there are a couple of potential problems.  I will go through each one and explain how you can overcome them.

What do you do if your data is not a nice multiple of the decimation factor?  You can easily overcome this by padding the end of your data with enough zeros to make it a multiple of the decimation factor.  The data is padded BEFORE it is FFT'ed.
Even though the method I explained is very simple, it is also non-ideal in that it can introduce ringing and other nasty artifacts in the time domain.  You can avoid that by filtering the frequency domain data before dropping the high frequency data.  You do this by FFT'ing your filter of length $l$, padding your data (before FFT'ing it) with at least $l-1$ zeros (please note that the number of data samples and the number of padding samples must BOTH be a positive multiple of the decimation factor- you can increase the padding length to meet this constraint), FFT'ing the padded data, multiplying the frequency domain data and filter, and then aliasing the high frequency (> 8 kHz) results down into the low frequency (< 8 kHz) results before dropping the high frequency results.  Unfortunately, since filtering in the frequency domain is a big topic in its own right, I won't be able to go into more detail in this answer.  I will say, though, that if you filter and are processing the data in more than one chunk, you will need to implement Overlap-and-Add or Overlap-and-Save to make the filtering continuous.

I hope this helps.

EDIT:  The difference between the starting number of frequency domain samples and the target number of frequency domain samples needs to be even so that you can remove the same number of samples from the positive side of the results as the negative side of the results.  In the case of our example, the starting number of samples was the decimation rate, or 441, and the target number of samples was the interpolation rate, or 160.  The difference is 279, which is not even.  The solution is to double the FFT length to 882, which causes the target number of samples to also double to 320.  Now the difference is even, and you can drop the appropriate frequency domain samples with no problems.


Do JPEG formats add significant overhead?

The picture you've posted most certainly didn't come directly from the image sensor – it was processed by the camera hardware, and software; it's pure dark grey.

Anyway; JPEG (to be more specific: JFIF as the JPEG-carrying file format) is not zero overhead, and you can expect a couple kB for this kind of image to work.

There's a couple bytes of headers (it's really negligible), but then things get a bit more complicated:

For a monochrome image, the actual JPEG quality setting won't do anything to the file size; it just sets a quantization mask, and you're always hitting the DC bin with a value of nonzero for anything. No matter how that turns out, Run-Length Encoding for the rest of the bins in each 8x8 pixel block will reduce the rest to the same number of bytes, regardless of quantization.
What does matter is the subsampling mask for the channels. Stanley's answer is a tiny bit wrong w.r.t. to the three channels: JPEG does support 4:2:0 subsampling, meaning that you actually omit one channel every row. But that's a story for a different day. Anyways, since your image doesn't change spatially at all, that saves space
Entropy encoding (I think it's a Huffmann method) is applied to the rest. Since the picture is built from identical 8x8 blocks, these should be reduced  to very little data at all.
There's the special JPEG/JFIF feature of interleaving the info in a way that allows the lower frequency parts of the image to be transferred first (e.g. for slow connections; if you grew up in the 90s, you'll remember pictures that got sharper gradually while you hoped for your dialup modem to take less costly time). While these practically make little difference for "typical" photos, for your monochrome picture, they'll make matters worse and increase the file size by some 20 to 25%, maybe.

So, 6 kB is actually pretty OK, and you should probably look into whatever gave you that picture – it's definitely heavily smoothing the image, to a degree that the data you're presented has nothing to do with what the sensor sees. The sensor sees noise.


how can I design Chebyshev filter with 5th-order polynomial and 0.001%-allowed ripples?

A type I Chebyshev filter has a passband magnitude response oscillating between the values $1$ and $1-\delta$, where $\delta$ is the maximum passband approximation error. The Matlab function $\tt{cheby1}$ requires the desired passband ripple in dB:

$$R_p=20\log\left(\frac{1}{1-\delta}\right)$$

The normalized cut-off frequency is

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

where $f_c$ is the cut-off frequency in Hertz, and $f_s$ is the sampling frequency in Hertz. The input parameter $\tt{Wc}$ required by $\tt{cheby1}$ is equal to $\omega_c/\pi$.

I don't think that the authors of the paper you refer to got the percentages right. A tip: for a "ripple percentage" of $0.001$ try $\delta=0.0001$.


Matlab - Segment cell nucleus from similar background

You need to work with different color space. Try YCbCr for instance

oimg = imread('http://i.stack.imgur.com/aYhrS.png');
img = rgb2ycbcr(oimg); % conver RGB to YCbCr color space
msk = img(:,:,3)>130;  % apply threshold on Cr channel to get segmentation mask

I picked up 130 as arbitrary threshold, you might need to adjust it a bit.
Showing the input image, the image region for which msk==1, and the background:

I'll leave it to you as an exercise to think why the Cr channel performed well for this image.


What is the total impulse response in a system with feedback interconnection?

It is easier to work in the $s$-domain:

$$Z=YH_2$$
$$Y=(X-Z)H_1$$
Hence,
$$Y=(X-YH_2)H_1=XH_1-YH_1H_2\Rightarrow Y(1+H_1H_2)=XH_1$$
Therefore,
$$H(s)=\frac{Y(s)}{X(s)}=\frac{H_1(s)}{1+H_1(s)H_2(s)}$$
which is called the closed-loop transfer function. The closed-loop impulse response can be found by inverse Laplace transform.


Why we sometime use some values of fixed point and floating point both in dsp?

The transition from text-book math to product code can be a lot simpler with floating point than with fixed point. The latter requires some specialized knowledge that, while not rocket science, may be more common for electronic engineers or embedded engineers than computer science or math people. Simple stuff tends to be a lot simpler in floating point than in fixed point.
If you really need to know the numerical behaviour of your code (as in certain recursive algorithms), it tends to be the other way around - the «cruise control» of floating point also makes it more tedious to figure out what is actually going on, while the explicit requirements of fixed point makes the details more plain sight.


Filtering out a specific frequency in an analog signal

A notch filter is a good approach in general for removing an unwanted frequency. But if you are certain of the exact frequency, phase and amplitude of the unwanted sine wave, then you may also consider adding the noisy signal to a sine wave of the same frequency and amplitude but 180 degrees out of phase with the unwanted sine wave.

To illustrate the notch filter example, the plot below shows an attempt to generate a sine wave with some odd impulses (black), and then filter it with a notch filter knowing the exact frequency (red).

R code for implementing it below.

# 24782
# install.packages('signal')

T <- 1024
t <- 0:(T-1)

fs <- 1000
omega <- 2*pi*30/1000
phi <- 0.9279835

x <- sin(omega*t + phi)

changed <- sample(1:T, 10)

y <- x
y[changed] <- y[changed] + 3

alpha <- 0.9
num <- c(1, -2*cos(omega), 1)
den <- c(1, -2*alpha*cos(omega), alpha*alpha)

yf <- signal::filter(num, den, y)

plot(t,y, type="l", lwd=10)
lines(t,yf, col="red", lwd=4)

legend(850, 4, c("Original", "Filtered"), lwd=c(2.5,2.5),col=c("black","red"))


Extending Goertzel algorithm to 24 kHz, 32 kHz and 48 kHz

Choice of $N$

Dual-tone multi-frequency (DTMF) tones should be of certain minimum duration (40 ms) and quality to be detected (reviewed in ITU TABLE A-1/Q.24). The detection algorithm that you cite uses a block duration of $t=N/f_s,$ where $N$ is the block length in samples and $f_s$ is the sampling frequency. The duration is either 210/(16 kHz) = 13.125 ms or 92/(8 kHz) = 11.5 ms. To ensure that the algorithm functions as expected, the block duration should be kept the same. I don't know why the lengths differ for the different sampling frequencies in the cited algo, but it should be okay to calculate $N$ as:

$$N = t\times f_s = 11.5\text{ ms}\times f_s,\tag{1}$$

which should give the same timing and frequency resolution as with the 8 kHz sampling frequency used in the original sources on Wikipedia by an unknown author.

Other changes needed

There are some worrisome large hard-coded constants (4.0e5 and 1.0e9) in the source code. Such hint to a lack of normalization:

#Check for minimum energy
if self.r[row] < 4.0e5:
    msg = "energy not enough"
elif self.r[col] < 4.0e5:
    msg = "energy not enough"
else:

#AT&T states that the noise must be 16dB down from the signal.
# Here we count the number of signals above the threshold and
#there ought to be only two.
if self.r[max_index] > 1.0e9:
    t = self.r[max_index] * 0.158
else:
    t = self.r[max_index] * 0.010

They may require changing if you change $N,$ by the same proportion, I think.

Fine-tuning $N$

Even if the Goertzel algorithm is not restricted to frequencies of a discrete Fourier transform (DFT) of length $N,$ (see this answer) there may be reason for tuning $N$ that is not explained by duration specifications. For sinusoidal input of frequency $f$, the magnitude of the same frequency detected by a Goertzel algorithm fluctuates slightly as function of time depending on the combination of frequency $f$ and $N$. For simplicity of the expressions that follow, frequency is expressed in radians as:

$$\omega=\frac{2\pi f}{f_s}\tag{2}$$

For initial phase $\alpha$ at the start of the window, the unnormalized detected magnitude will be:

$$A(\alpha) = \left|\sum_{k=0}^{N-1} e^{i \omega k}\cos(\omega k + \alpha)\right|\\
= \left|\frac{Ne^{-i\alpha}}{2} + \frac{e^{i(\alpha - \omega+\pi/2)} - e^{i(\alpha + 2N\omega - \omega + \pi/2)}}{4\sin(\omega)}\right|\tag{3}$$

The fluctuation can be seen in this example:

Figure 1. Blue: detected magnitude $A$ as function of initial phase $\alpha$ for sinusoidal input, with $f = 1209\text{ Hz},$ $f_s = 8 kHz,$ $N=92$. Red: extrema of $A(\alpha)$ have been marked.

An interesting observation is that there will no fluctuation of the Goertzel detected magnitude $A$ if the frequency $\omega$ is a multiple of $\pi/N.$ In contrast, DFT bins are at multiples of $2\pi/N.$ (Zaplata & Kasal 2005 discuss the fluctuation phenomenon, but do not recognize this difference between the Goertzel algorithm and DFT.)

For other frequencies, the $\pi\text{-}$periodic function $A(\alpha)$ gets its smallest and largest values (not necessarily in that order) at:

$$\alpha_1 = \frac{\omega(1 - N)}{2}\quad\text{and}\quad\alpha_2 = \frac{\omega(1 - N)}{2} + \frac{\pi}{2}\tag{4}$$

You could tune the integer $N$ to its nearby value which minimizes the maximum of $\left|A(\alpha_1)-A(\alpha_2)\right|$ over the set of DTMF frequencies. Actually, for $f_s = 8\text{ kHz,}$ a better choice than $N=92$ would be $N=93:$

Figure 2. Maximum peak-to-peak fluctuation $\max\left|A(\alpha_1)-A(\alpha_2)\right|$ over DTMF frequencies has a minimum at $N=93$ for $f_s=8000\text{ Hz.}$

A reason why $N=92$ might have been preferred is that in some older version they only allowed DFT bin frequencies. In that case $N=93$ would be a pretty bad choice, as some DTMF frequencies would fall almost just between successive bins, which is perfectly fine for the Goertzel algorithm (or generalized DFT with a frequency shift of $1/2$ binwidth), but almost the worst-case scenario when using DFT:

Figure 3. Maximum absolute frequency error over DTMF frequencies when detecting using DFT of length $N$, for $f_s=8000\text{ Hz}$ and expressed in units of DFT bin width.

The final word would be to test the performance in a test suite and to tune $N$ based on that.

Octave source code for Fig. 3

graphics_toolkit("gnuplot");
fs = 8000; #Sampling frequency
f = [1209, 1336, 1477, 1633, 697, 770, 852, 941]; #DTMF frequencies (same unit as fs)
N = 80:100; #DFT lengths analyzed
omega = 2*pi*f/fs;
y=omega'/(2*pi).*N;
z=round(y)-y;
plot(N, max(abs(z)), "x");
ylim([0, 0.5]);
xlabel("N");
ylabel('max|frequency error| (bin widths)');

References

Zaplata F, Kasal M, "Efficient Spectral Power Estimation on an Arbitrary Frequency Scale", Radioengineering, 2015, vol. 24: p. 178-184. DOI: 10.13164/re.2015.0178


Fourier transform of Image to identify sinusoidal sources of interference

Alright. This is image processing. Sometimes its results might seem indeed magical. And quite some other times, in some books or papers, you see they discuss an algorithm or method to process a given distorted/noisy image and convert it into an unbelievably perfected form. Yet they don't provide any actual code, neither are the input and output images made available. It's therefore quite impossible to replicate what they propose to have achieved.

In this example a ridiculously simple idea is provoked. In fact it is the most basic idea of all signal procesing: Filtering! There is an original degraded image so much degraded that it might seem impossible to recover, but when it's analysed via DFT it's been seen that there're some spurious frequency components which are not to be expected from the original undistorted image. In addition it's quite interesting that those spurious frequencies are all distributed inside a circular ring of frequencies! So this is the easiest of all signal processing. Just design a circular 2D filter narrow enough to suppress all those circular band of spurious components, and leave all the rest untuoched. The result is that magic.

That being said, however, It is very unlikely that in your everyday life you would encounter such a highly structured "noise" or "distortion". It seems as if almost this distortion was just added there for the purpose of a filtering demonstration. And indeed due to poor fact checking and documentation this kind of demonstration tricks seem to exist. What you basicly do is this: You have an original clear image. You add some artificially generated noise to it and you apply a noise recovery algorithm to remove it! Aaaand ta daaaa perfect recovery is possible, in real life or just for an artifically generated noise?


Personally I find Python one of the best choices out there and did myself some work in area of audio identification. You are welcomed to check for instance my software for automatic identification of birds from noisy audio recordings: Ornithokrites. The program is used by Department of Conservation of New Zealand and they are happy about it. Based on this example I would like to point out several advantages of using Python:

Huge, fast developing community providing tons of libraries. SciPy provides plethora of methods for signal processing (granted, not that many and mature as Matlab). Mind though that SciPy, although one of the most important, is only one of hundreds that can help you in your endeavours. I found Aubio best for music analysis. For speech and music recognition for sure you will enjoy great number of audio features Yaafe can extract.
It's free! Once out of academia, you quickly find out that Matlab is rather expensive. And even if you can afford it, then your perspective users will not be happy about this dependency. For instance, mentioned Department of Conservation would not accept proprietary software.
Identification often requires machine learning and Python has great toolkit for it: sklearn. It is state of the art library - and easy to use. Have a look at Kaggle competitions  (machine learning) and check how many top programmers are using Python and sklearn.
You can manage "big data". If you want to run analysis against huge networked databases of recordings, then Python has well established set of tools. I don't think Matlab / Octave  interface easily with e.g. Hadoop, although please do correct me if I am wrong. R does better this area.
Speaking of interfacing, you can easily interface your program with a web site. This is the way I manage Ornithokrites (bird recognition): the program runs on Amazon Web Services cloud computing service. Great if you want to provide your software to other people who do not necessarily want to go through installation procedure of all required libraries.

My second choice would be R. Although not that feature-rich as Python, it has great number of useful libraries (check e.g. seewave for your applications). Installation of those on both Windows and Linux is piece of cake, which is important if you would like others to use your program. However, to my experience high-performance computing in R is more difficult - an important thing to notice if you need to do A LOT of processing and identification.

Examples of music classification in Python:

Lecture
Audio fingerprinting
Jack Minardi blog

Book Building Machine Learning Systems with Python has a chapter on music classification

Other tools (list by no means complete): Python in Music


What is an algorithm to re-sample from a variable rate to a fixed rate?

The simplest approach is to do some kind of spline interpolation like Jim Clay suggests (linear or otherwise). However, if you have the luxury of batch processing, and especially if you have an overdetermined set of nonuniform samples, there's a "perfect reconstruction" algorithm that's extremely elegant. For numerical reasons, it may not be practical in all cases, but it's at least worth knowing about conceptually. I first read about it in this paper.

The trick is to consider your set of nonuniform samples as having already been reconstructed from uniform samples through sinc interpolation. Following the notation in the paper:

$$
y(t) = \sum_{k=1}^{N}{y(kT)\frac{\sin(\pi(t - kT)/T)}{\pi(t - kT)/T}} = \sum_{k=1}^{N}{y(kT)\mathrm{sinc}(\frac{t - kT}{T})}.
$$

Note that this provides a set of linear equations, one for each nonuniform sample $y(t)$, where the unknowns are the equally-spaced samples $y(kT)$, like so:

$$
\begin{bmatrix} y(t_0) \\ y(t_1) \\ \cdots \\ y(t_m) \end{bmatrix} = \begin{bmatrix} \mathrm{sinc}(\frac{t_0 - T}{T}) & \mathrm{sinc}(\frac{t_0 - 2T}{T}) & \cdots &  \mathrm{sinc}(\frac{t_0 - nT}{T}) \\ \mathrm{sinc}(\frac{t_1 - T}{T}) & \mathrm{sinc}(\frac{t_1 - 2T}{T}) & \cdots & \mathrm{sinc}(\frac{t_1 - nT}{T}) \\ \cdots & \cdots & \cdots &\cdots \\ \mathrm{sinc}(\frac{t_m - T}{T}) & \mathrm{sinc}(\frac{t_m - 2T}{T}) & \cdots & \mathrm{sinc}(\frac{t_m - nT}{T}) \end{bmatrix} \begin{bmatrix} y(T) \\ y(2T) \\ \cdots \\ y(nT) \end{bmatrix}.
$$

In the above equation, $n$ is the number of unknown uniform samples, $T$ is the inverse of the uniform sample rate, and $m$ is the number of nonuniform samples (which may be greater than $n$).  By computing the least squares solution of that system, the uniform samples can be reconstructed. Technically, only $n$ nonuniform samples are necessary, but depending on how "scattered" they are in time, the interpolation matrix may be horribly ill-conditioned. When that's the case, using more nonuniform samples usually helps.

As a toy example, here's a comparison (using numpy) between the above method and cubic spline interpolation on a mildly jittered grid:

(Code to reproduce the above plot is included at the end of this answer)

All that being said, for high-quality, robust methods, starting with something in one of the following papers would probably be more appropriate:

  A. Aldroubi and Karlheinz Grochenig, Nonuniform sampling and
  reconstruction in shift-invariant spaces, SIAM Rev., 2001, no. 4,
  585-620. (pdf link).

  K. Grochenig and H. Schwab, Fast local reconstruction methods for
  nonuniform sampling in shift-invariant spaces, SIAM J. Matrix Anal.
  Appl., 24(2003), 899-
  913.

--

import numpy as np
import pylab as py

import scipy.interpolate as spi
import numpy.random as npr
import numpy.linalg as npl

npr.seed(0)

class Signal(object):

    def __init__(self, x, y):
        self.x = x
        self.y = y

    def plot(self, title):
        self._plot(title)
        py.plot(self.x, self.y ,'bo-')
        py.ylim([-1.8,1.8])
        py.plot(hires.x,hires.y, 'k-', alpha=.5)

    def _plot(self, title):
        py.grid()
        py.title(title)
        py.xlim([0.0,1.0])

    def sinc_resample(self, xnew):
        m,n = (len(self.x), len(xnew))
        T = 1./n
        A = np.zeros((m,n))

        for i in range(0,m):
            A[i,:] = np.sinc((self.x[i] - xnew)/T)

        return Signal(xnew, npl.lstsq(A,self.y)[0])

    def spline_resample(self, xnew):
        s = spi.splrep(self.x, self.y)
        return Signal(xnew, spi.splev(xnew, s))

class Error(Signal):

    def __init__(self, a, b):
        self.x = a.x
        self.y = np.abs(a.y - b.y)

    def plot(self, title):
        self._plot(title)
        py.plot(self.x, self.y, 'bo-')
        py.ylim([0.0,.5])

def grid(n): return np.linspace(0.0,1.0,n)
def sample(f, x): return Signal(x, f(x))

def random_offsets(n, amt=.5):
    return (amt/n) * (npr.random(n) - .5)

def jittered_grid(n, amt=.5):
    return np.sort(grid(n) + random_offsets(n,amt))

def f(x):
    t = np.pi * 2.0 * x
    return np.sin(t) + .5 * np.sin(14.0*t)

n = 30
m = n + 1

# Signals
even   = sample(f, np.r_[1:n+1] / float(n))
uneven = sample(f, jittered_grid(m))
hires  = sample(f, grid(10*n))

sinc   = uneven.sinc_resample(even.x)
spline = uneven.spline_resample(even.x)

sinc_err   = Error(sinc, even)
spline_err = Error(spline, even)

# Plot Labels
sn = lambda x,n: "%sly Sampled (%s points)" % (x,n)
r  = lambda x: "%s Reconstruction" % x
re = lambda x: "%s Error" % r(x)

plots = [
    [even,       sn("Even", n)],
    [uneven,     sn("Uneven", m)],
    [sinc,       r("Sinc")],
    [sinc_err,   re("Sinc")],
    [spline,     r("Cubic Spline")],
    [spline_err, re("Cubic Spline")]
]

for i in range(0,len(plots)):
    py.subplot(3, 2, i+1)
    p = plots[i]
    p[0].plot(p[1])

py.show()


why z-transform differentiation is needed?

You just wrote down the $\mathcal{Z}$-transform of $\left(\frac12\right)^nu[n]$, but you need the $\mathcal{Z}$-transform of $n\left(\frac12\right)^nu[n]$ (note the factor $n$).

In order to find that $\mathcal{Z}$-transform you can use the differentiation property (see this table):

$$\mathcal{Z}\big\{nx[n]\big\}=-z\frac{dX(z)}{dz}\tag{1}$$


Uniform quantizer for gaussian input signal

The problem with quantizing Gaussian distributed signals (like the real/imaginary part of an OFDM signal) is that they can take any value in theory. It is thus necessary to clip such signals at threshold $C$ prior to quantization.

Low $C$ increases the distortion noise in this process, while large $C$ will lead to strong quantization noise. It can be shown, that there is an optimum value for $C$ that depends on the variance $\sigma^2$ of the signal and the quantizer resolution $M$. For details see [1,2]. After $C$ has been determined, the step size $q$ is calculated by
$$
q=\frac{2C}{2^M}
$$

By the way, I don't think it makes sense to use the ENOB (of a DAC) for calculating the step size. The ENOB is a kind of signal-to-noise ratio combining quantization noise and other effects like frequency selectivity and spurious frequencies of the DAC device. The number of input bits is still the nominal number of bits $M$.

Example: let $\sigma^2=1$ and $M=8$. Then with eq. (7) from [2] we can calculate the SNR $\gamma$ after clipping and quantization as a function of the clipping ratio $\zeta=C^2/\sigma^2$. The SNR should be maximized. I don't know a closed-form expression for the optimum $\zeta_\mathrm{opt}$ yielding $\gamma_\mathrm{max}$ but we can determine it numerically. For $M=8$ and with Fig. 2 from [2] we have $\zeta_\mathrm{opt} \approx 11.87\,\text{dB} \approx 15.38$. Furthermore, $C_\mathrm{opt}=\sqrt{\zeta_\mathrm{opt}\sigma^2}$ and with $\sigma^2=1$ we finally get $C_\mathrm{opt}\approx3.92$.

Remark: as you tagged this question digital-communictions you might be dealing with complex signals. In that case be careful with the definition of $\sigma^2$. In my answer it is the variance of a real-valued, Gaussian distributed signal. If the signal in question is complex-valued with independent and identically Gaussian distributed real and imaginary part, it has variance $2\sigma^2$, where $\sigma^2$ is the variance of the real or imaginary part.

References

[1] D. J. G. Mestdagh, P. Spruyt, and B. Biran, “Analysis of clipping effect in DMT-based ADSL systems,” in Communications, 1994. ICC’94, SUPERCOMM/ICC'94, Conference Record,'Serving Humanity Through Communications.'IEEE International Conference on, 1994, pp. 293–300.

[2] M. Bernhard, D. Rörich, T. Handte, and J. Speidel, “Analytical and numerical studies of quantization effects in coherent optical OFDM transmission with 100 Gbit/s and beyond,” in ITG Symposium on Photonic Networks, 2012, pp. 34–40. Online


Let's have a different perspective on that.

Let's say our 2D Linear Operator is given by the Matrix $ G \in {\mathbb{R}}^{n \times n} $.
Using the SVD Decomposition the operator can be written as:

$$ G = \sum_{i = 1}^{n} {\sigma}_{i} {u}_{i} {v}_{i}^{T} $$

Separable Linear 2D Operator is defined as operator which can be composed by Outer Product of 2 vectors.
Looking at the SVD Decomposition of $ G $ we can conclude that $ G $ is separable operator if and only if $ \forall i > 1 \; {\sigma}_{i} = 0 $ and it is given by:

$$ G = {\sigma}_{1} {u}_{1} {v}_{1}^{T} $$

Usually LPF 2D Linear Operators, such as the Gaussian Filter, in the Image Processing world are normalized to have sum of 1 (Keep DC) which suggests $ {\sigma}_{1} = 1 $ moreover, they are also symmetric and hence $ {u}_{1} = {v}_{1} $ (If you want, in those cases, it means you can use the Eigen Value Decomposition instead of the SVD).

So basically, to prove that a Linear 2D Operator is Separable you must show that it has only 1 non vanishing singular value.

Note: I was not rigorous in the claims moving form general SVD to the Eigen Decomposition yet the intuition holds for most 2D LPF operators in the Image Processing world.


How can you tell whether there's a signal (when your signal looks a lot like noise)?

So far you seemed to identify snoring through the presence of a strong periodic pattern in the audio. If you tell me that there might be other sources with the same property, it's time to move on and focus on a property of the signal more specific of snoring ; and I would suggest to look more precisely into the timbre of the recorded sound itself. Given that a human listener can tall apart snoring from breathing, you might resort to traditional audio signal classification approach. Compute features like MFCCs and Fluctuation patterns on your audio data, and train a bayesian classifier (or more sophisticated) to discriminate snoring from breathing.

It might be time also to review the literature on the topic... For example I found this: http://web.itu.edu.tr/kamasak/pubs/pdf/pm7_8_007.pdf


Why is dirac delta used in continuous signal sampling?

The best explanation of this that I've seen is from the "Digital Signal Processing Handbook" by Madisetti. Essentially the multiplication by the delta function is equivalent to sampling because the Fourier transforms are the same. So although the result of $s(t)\delta(t-nT)$ may not make much sense it's Fourier transform does exist because we are taking the integral over the delta function.

I'm cutting and pasting his text and redoing the equations. Note CT=continuous time, DT= discrete time. I have not included figures 1.1 and 1.2.

The relationship between the CT and the DT domains is characterized by the operations of sampling and reconstruction. If $s_a(t)$ denotes a signal $s(t)$ that has been uniformly sampled every T seconds, then the mathematical representation of $s_a(t)$ is given by

$s_a(t) =\sum_{n=- \infty}^{\infty} s(t)\delta(t-nT) \qquad$   (1.1)

where $\delta(t)$ is a CT impulse function defined to be zero for all $t\not=0$, undefined at t=0, and has
unit area when integrated fromt $t=-\infty$ to $t=\infty$. Because the only places at which the product $s(t)\delta(t−nT)$ is not identically equal to zero are at the sampling instances, $s(t)$ in (1.1) can be replaced with $s(nT)$ without changing the overall meaning of the expression. Hence, an alternate expression for $s_a(t)$ that is often useful in Fourier analysis is given by

$s_a(t) =\sum_{n=- \infty}^{\infty} s(nT)\delta(t-nT) \qquad$       (1.2)

The CT sampling model $s_a(t)$ consists of a sequence of CT impulse functions uniformly spaced at intervals of $T$ seconds and weighted by the values of the signal $s(t)$ at the sampling instants, as depicted in Fig.1.1. Note that $s_a(t)$ is not defined at the sampling instants because the CT impulse function itself is not defined at $t$. However, the values of $s(t)$ at the sampling instants are imbedded as “area under the curve” of $s_a(t)$, and as such represent a useful mathematical model of the sampling process. In the DT domain the sampling model is simply the sequence defined by taking the values
of $s(t)$  at the sampling instants, i.e.,

$s[n]=s(t)|_{t=nT}  \qquad $(1.3)

In contrast to $s_a(t)$, which is not defined at the sampling instants, $s[n]$ is well defined at the sampling instants, as illustrated in Fig.1.2. Thus, it is now clear that $s_a(t)$ and $s[n]$ are different but equivalent models of the sampling process in the CT and DT domains, respectively. They are both useful for signal analysis in their corresponding domains. Their equivalence is established by the fact that they
have equal spectra in the Fourier domain, and that the underlying CT signal from which $s_a(t)$ and $s[n]$ are derived can be recovered from either sampling representation, provided a sufficiently large sampling rate is used in the sampling operation (see below).

Cheers,

David


you can use this approach provided it does not yield efficiency burdens.

% Note that matlab arrays start from index 1, not 0.
F1 = 223807.5;     % Sampling rate of Yamaha chip
F2 = 44100;        % New sampling rate of standard Audio-CD
M = 203;           % Necessary Decimation ratio
N = 40;            % necessary Expansion ratio
Lx = 4466;         % Length of Audio Input Buffer at F1 rate
Ly = Lx*N/M;       % Length of Resampled Output Buffer at F2 rate

x = sin(2*pi*1234*[0:Lx-1]/F1);  % Just a test sine wave
y = zeros(1,Ly);                 % allocate output signal memory

for i=0:Ly-1              % run per each output sample
   n = i*(M/N);          % compute exact sampling position inside the long buffer
   ind = floor(n);       % convert it into largest integer smaller than itself for array indexing
   d = n-ind;            % take the difference for proportional weighting
   y(i+1) = (1-d)*x(ind+1) + d*x(ind+2); % sum the weighted mixture.
end

figure,plot(x);            %plot each waveforms...
figure,plot(y);


This is what I do when making such a plot for a QPSK modulation:

First I upsample the I and Q data so that I can see more of the trajectory between samples given (assuming you start with at least 2 samples per symbol otherwise, due to Nyquist, the information on the trajectory between symbols is not available, but as long as you have 2 or more samples, all the information of what happens in the first Nyquist zone is preserved):

Is = resample(I, 10, 1);

Qs = resample(Q, 10, 1);

and then I make a plot of I vs Q

plot(Is,Qs)

I used a resampling increase of 10, but you could make it higher if desired. The result is the interpolated samples for the full trajectory during symbol transitions.


Harry Nyquist made so many contributions that it's easy to get confused.

Related to sampling, Nyquist proved that a signal $s(t)$ bandlimited to $B$ Hz can be reconstructed from samples taken at a rate larger than $2B$. This is unrelated to Faster than Nyquist (FTN) signaling, though.

Related to communications, Nyquist also showed that the maximum  ISI-free pulse rate in a channel of bandwidth $B$ is $2B$ pulses per second, using a technique known as "orthogonal signaling".

Note that Nyquist didn't say that communication at a rate faster than $2B$ was impossible; he just said that you will have to handle the ISI somehow. For a long time, people didn't know how to do that, so $2B$ became engrained as a sort of actual bound.

One exception is partial response signaling, which introduces ISI in a controlled manner. However, in this case the main purpose is spectral shaping more than an increase in data rate.

In the 1970s, Mazo published a paper showing that, in some cases, the ISI can not only be handled, but it actually doesn't detract from a system's error rate performance. Start from an orthogonal signaling system operating at pulse rate $R=1/T$; that is, find an orthonormal set of pulses $p(t-kT)$ (for integers $k,l,m$) such that
$$
p(t-lT)p(t-mT) = \begin{cases}
1,\text{when $l=m$} \\
0,\text{otherwise.}
\end{cases}
$$
Now let's say that you want to transmit the data sequence $\lbrace a_0, a_1, a_2, \ldots, a_M \rbrace$. You can do that by transmitting
$$
s(t) = \sum_{k=0}^M a_k p(t-kT)
$$
and the receiver (ignoring noise for simplicity) can recover any $a_k$ by calculating
$$
a_k = \int_{-\infty}^{\infty} s(t) p(t-kT)\, dt.
$$

Mazo's proposal was to convert this orthogonal-signaling system to an FTN system by introducing a time-compression factor $0 < \tau < 1$, so that the transmitted signal becomes
$$
s(t) = \sum_{k=0}^M a_k p(t-k \tau T)
$$
Note that the bandwidth of $s(t)$ is still $B$, but now we transmit pulses at rate $1/(\tau T) > R$. Mazo's main result was to prove that, for certain pulses $p(t)$ and certain time-compression factors $\tau$, the bit-error rate performance of the accelerated system is the same as the orthogonal system.

At the time, the problem was that handling the ISI was impractical. The ISI acquires a certain structure that can be represented as a trellis, but in any non-trivial case the trellises were huge and well beyond the computational capabilities available.

It is only recently that some new techniques, along with faster computers, have made FTN somewhat practical. I recommend this tutorial paper:

J. B. Anderson, F. Rusek, and V. wall, “Faster-Than- Nyquist Signaling,” Proceedings of the IEEE, vol. 101, no. 8, pp. 1817–1830, Aug. 2013.

Mazo's original paper is:

J. E. Mazo, “Faster-Than-Nyquist Signaling,” Bell System Technical Journal, vol. 54, no. 8, 1975.


RMS includes the mean. If you want a metric for "variation" use the standard deviation instead or subtract the mean first (which is basically the same thing).
For the black curves moves slowly, but it still moves by a good amount. If you want a metric of the "fuzziness" than you can apply a high pass filter first that gets rid of the slow "drift" movemenent.


3D (time, scale, amplitude) plot in Continuous Wavelet Transform

It suffices to use the translation and scale as $(X,Y)$ axes, and build some elevation map from the absolute values of the CWT (or the phase, the real or imaginary parts). A simple example in Matlab is:
load mtlb
scal = abs(cwt(mtlb,'bump',Fs));
mesh(scal)
axis tight


How to select the sampling frequency when the input signal frequency is not known

When the input signal bandwidth in unknown, you can begin with the highest possible sampling frequency, then observe the energy in the high frequency region (close to Nyquist frequency; i.e. half the sampling rate) and try to estimate the bandwidth of the signal and then choose the lowest possible rate still maintainig the Nyquist criterion according to this accepted bandwidth.

If the first trial indicates that there's considerable energy towards the Nyquist frequency, then you should try more complex sampling strategies, such as bandpass sampling at higher frequency regions, unless you can afford an even higher frequency sampling device.


How to achieve full number of words in a Bag-of-Features model based on Hierarchical K-Means?

Unfortunately, I think what causes the problems is most probably the training set of images you're using to construct the tree. (or, more in general, the database you're working on). Without a way to inspect the concrete images, I can only give a few vague guesses.

First, I would diagnose where it goes wrong, by simply inspecting the tree and looking at the first branches which become severely underpopulated and consequently don't branch any more in the lower levels. I would look at those features and... think. But, since I can't actually do that since I don't have your dataset, let me try further.

If you're afraid that it's the K-means that is going wrong, not because the features aren't (more-or-less) uniformly distributed over the feature space, but simply because of the way that algorithm works, before implementing a fancier algorithm, you might first try running K-means several times for each branching and picking the best result (e.g. the one that divides the feature space in most equal parts). This will probably be too slow for the final implementation, but it could be a fast, minimum-effort way do diagnose the problem.

I would say, in order for HKM to work, your initial number of features should be several orders of magnitude greater than the expected number of clusters. So, if that's not the case, you know you're doing something wrong. I don't know how many is enough, but I would say minimum $100$ times, but better at least $1000$ times more features than expected clusters. This might be caused by two (related) issues:

Maybe your training set is just too small, or maybe it is a bad representation of the database, containing images that are too similar? Try increasing the size of the training set, and maybe (just for testing purposes), cherry pick the training set by selecting very different, representative images.
Maybe your feature detector performs badly on the image set you're using. For some images, some detectors just do not detect enough features. Problematic images usually contain large single-colored areas (e.g. white wall, blue sky). Try combining several feature detectors. (e.g. in my thesis I was working with DoG and MSER features). You can always do that, since you can finally describe all the features from several detectors with the same descriptor.

Unfortunately, there's no more ideas I can think of based on so little information. You might simply have an implementation error sometimes. I know it's not easy to here, and it's not a most helpful diagnosis, but as somebody who actually implemented such a CBIR system, I have to tell you, I never had any problem in constructing a vocabulary tree of expected size.

If you give any more info, or updates, or whatever, I can try and think of other ways this might have gone wrong. Until then, I hope this helps.


Least Squares Filter Design: Deriving the Objective Function

$$
J(h) = \int_{R^+UR-}|e^Th-F_d|^2d\omega\tag{1}\\
= \int(e^Th-F_d)^H(e^Th-F_d)d\omega\\
= \int((e^Th)^H(e^Th) + F_d^HF_d -(e^Th)^HF_d - F_d^He^Th)d\omega\\
= h^H(\int(e^T)^He^Td\omega) h + \int |F_d|^2d\omega -\int (2 Re\{(e^Th)^HF_d\})d\omega
$$
The second term in above integral is the integral of $L_2$ norm of $F_d$ in the region $R^+ U R^-$. Since $F_d(\omega)$ is symmetric about $\omega = 0$, this integral will be $2 \times \int_{R^+}|F_d(\omega)|^2d\omega$ which is $c$ in your questions.

The first term, the argument inside integral is the outer product of column vector $(e^H)^T$ with row vector $e^T$. Again, this will conjugate of each for the regions $R^+$ and $R^-$ (shown later in the appendix). So integral over both regions will result in the real value of integral $e^H e$.  Also, since $H$ is real, $h^H=h^T$. So the first term can be rewritten as (for simplicity dropping $T$ notation as it is assumed to be a column vector)
$$
h^T(\int_{R^+}2Re\{ ee^H\}d\omega )h\\
= h^TQh
$$

The third term is the result of sum of 2 conjugate terms. Again this is symmetric about $\omega=0$. Also, since $H$ is real, $h^H=h^T$. Taking out $2h^T$ from the integral we have
$$
2h^T\int_{R+}2Re\{F_de^H \}d\omega = 2h^Tp
$$
Therefore, summing and rearranging the 3 terms,
$$
J(h) = h^TQh - 2h^Tp + c
$$

Appendix: Showing that $ee^H$ is conjugate of each other in $R^+$ and $R^-$. $E_{mn}$, element of $ee^H$ is $e^{-jm\omega}e^{-jn\omega} = e^{-j(m+n)\omega}$. In $R^-$, this element is $e^{j(m+n)\omega}$ which is $E_{mn}^*$


Is there any research work on the possible data structures that can be used to represent an architectural design of a building?

DXF in computer aided design should already allow this representation. I would suggest you to check the format and adapt your data structures accordingly.

Check the data structures in LibDWG and LibDXF.


[1] It's a notch filter. It attenuates everything but the frequency component you're interested in. The output is summed up over some time to make the tone detection more reliable.
[2] There's no minimum block size. There's a requirement about the amount of time spent summing up the output. What you'll see is that the summed output over say 205 samples (at fs 8000kHz), grows much more if the input contains the frequency component you're interested in detecting.


Why can't a causal digital filter have an infinitely sharp transition between the passband and the stopband?

I don't have a concrete proof for this one. However, I can tell you this...
Consider a perfect low pass filter. The time domain representation is a sinc. And for any system to have a sharp transition band, a base signal has to be multiplied with a rectangular waveform in the frequency domain. Which implies that, the time domain signal of the same has to be convoluted with a sinc in time domain. We know that sinc function is a non causal signal.
Hence proved.


The EURASIP Journal on Advances in Signal Processing is a Springer journal with open access.  It sounds like the topic of your paper is appropriate for it.  The impact factor is not as high as IEEE Transactions on Signal Processing though.


What is the unit for power spectrum of acceleration signals in log scale

The typical practice with decibels is to append the unit.  Unfortunately, it's informal, so you'll see things like "dBm" in RF circuits which refers to "dB with a reference of 1mW", not "dB with a reference of 1 meter".

you could:

Put a note in your graph that $0\mathrm{dB} = \mathrm{1 m/s^2}$
Use "dBg" (and use 1g as your reference)
Just keep the y-axis as the log of acceleration squared (and yes, $\mathrm{m^2/s^4}$ is awkward, but it's the correct unit).


Design Lowpass Filter With -7.5 dB/Octave Rolloff

As far as I know there is no standard procedure for designing such filters. A simple approach would be to design a standard second-order low pass filter with the given cut-off frequency (such as a Butterworth filter), and then re-design the numerator coefficients such that the total filter approximates the desired magnitude roll-off. This is basically a (linear phase) FIR filter design problem, which is very straightforward to solve.

Let $M(\omega)$ be the desired magnitude response. On a log-log scale, $M(\omega)$ is piecewise linear ($0$ dB in the passband, linear decay $r$ dB/decade in the stopband). Furthermore, let $A(\omega)$ be the frequency response of the denominator coefficients of the original second-order filter. Then the desired response for the new numerator coefficients is given by

$$D(\omega)=M(\omega)\cdot |A(\omega)|\cdot e^{-j\omega}\tag{1}$$

where the exponential term is a linear phase term corresponding to the delay of a second-order linear phase FIR filter ($1$ sample).

Any method can be used to design that $3$-tap FIR filter (the numerator coefficients). I used a least squares design and the result is shown in the figure below (with the original Butterworth filter in blue, the desired magnitude response in red, and the modified Butterworth filter with re-designed numerator coefficients in green):

If you need a better approximation you can increase the order of the numerator (without changing the second-order denominator coefficients). So you'd get a biquad concatenated with a (short) FIR filter. The figure below shows the approximation with a $10^{th}$-order numerator polynomial and the second-order Butterworth denominator coefficients (in green):


This is solved slightly differently if you're simulating a continuous-time or discrete-time AWGN channel. I'll assume you're simulating a discrete-time channel; that is, you're simulating the output of the receiver's matched filter sampled at correct time instants.

First consider the case with no noise. In this scenario, the matched filter's output is an array of numbers, each of which is equal to $+\sqrt{E_b}$ or $-\sqrt{E_b}$. The average energy per sample is then equal to $E_b$.

Now consider the case where only white noise is present at the matched filter's input. The sampled filter output is a sequence of random values with variance $\sigma_n^2$. The signal-to-noise ratio is defined as $E_b/\sigma_n^2$. By convention, $\sigma_n^2$ is denoted $N_0/2$, which means that $\text{SNR}=2E_b/N_0$. (This is because when the filter input is white noise with power spectral density $N_0/2$ and the matched filter's impuse response has energy equal to 1, then the variance of the sampled filtered noise is also $N_0/2$. This is a technical detail that you may omit for now and come back to when you're more comfortable with the theory).

Then, $$\sigma_n^2=\frac{N_0}{2}=\frac{E_b}{SNR}.$$ Using this equation you can easily simulate any SNR you want. A word of caution: it is common, when plotting the bit error rate vs the SNR, to use $E_b/N_0$. Be aware of this when comparing your results to those published elsewhere.


Design causal FIR filter to approximate hilbert transform

Basically, Hilbert transformers are, by definition, non-causal, always. Any zero-phase filter is non-causal. We can shift the impulse in time but the phase then becomes linear and not flat (a function of frequency).
To be implemented, the input signal must be delayed by the amount that is half the length of the filter minus one (N in this case). Or one can convolve the input signal with the impulse of the Hilbert (windowed, as you did, to minimize Gibbs). Never done the latter so check yourself.
BUT! There is a configuration that uses a pair of 4th order allpass filters that turn out to be pretty close to 90 degrees out of phase with one another. Unfortunately, the phase of each output differs from the original but depending on what you need, it can be hacked...
Reference: https://www.dspguide.com/ch19/4.htm, Design Hilbert filter at matlab (pure 90 phase and magnitude 0db)


On computing the number of independent values from a signal's bandwidth?

It comes from the Nyquist rate.  The Nyquist sample rate for a signal with bandwidth $B$ is $2B$.  In other words, that is the lowest sample rate that will contain all of the information in the signal without distorting it (i.e. without losing information).

The question implicitly assumes that the sample rate is $2B$, and thus that in time $T$ you will get $2B*T$ independent samples.  If the sample rate was greater than $2B$, then the signal would not be "white" anymore and the samples would have some correlation.


You're right that the strictly recursive approach is vulnerable to the accumulation of error as the number of iterations increases. One more robust way this is typically done is to use a numerically-controlled oscillator (NCO). Basically, you have an accumulator that keeps track of the instantaneous phase of the oscillator, updated as follows:

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

$$
\phi[n] = (\phi[n-1] + \delta) \mod 2\pi
$$

At each time instant, then, you're left with converting the accumulated phase in the NCO to the desired sinusoidal outputs. How you do this depends on your requirements for computational complexity, accuracy, etc. One obvious way is to just calculate the outputs as

$$
x_c[n] = \cos(\phi[n])
$$

$$
x_s[n] = \sin(\phi[n])
$$

using whatever implementation of sine/cosine you have available. In high-throughput and/or embedded systems, the mapping from phase to sine/cosine values is often done via a lookup table. The size of the lookup table (i.e. the amount of quantization you do on the phase argument to sine and cosine) can be used as a tradeoff between memory consumption and approximation error. The nice thing is that the amount of computations required is typically independent of the table size. In addition, you can limit your LUT size if needed by taking advantage of the symmetry inherent in the cosine and sine functions; you only really need to store one-fourth of a period of the sampled sinusoid.

If you need higher accuracy than a reasonably-sized LUT can give you, then you can always look at interpolation between table samples (linear or cubic interpolation, for example).

Another benefit of this approach is that it is trivial to incorporate frequency or phase modulation with this structure. The frequency of the output can be modulated by varying $\delta$ accordingly, and phase modulation can be implemented by simply adding to $\phi[n]$ directly.


EDIT : Based on RBJ's comment that the indicated the system cannot possibly be an LTI, due to the fact that the output spectrum contains frequencies that do not exist at the input spectrum..., I've looked at the input-output spectrums, and came to the following conclusion :
Given the input:
$$ s(t) = \begin{cases} { 1 ~~~,~~~-0.5 < t < 0.5 \\ 0~~~,~~~\text{otherwise} }\end{cases} \tag{1}
$$
and the output:
$$ g(t) = \begin{cases} { 2t+1 ~~~,~~~ -0.5 < t < 0 \\ 1-2t ~~~,~~~ ~~~~~~~~0 < t < 0.5 }\end{cases} \tag{2}
$$
The corresponding spectrums (apart from linear scale factors !) are :
$$ S(\omega) = \frac{ \sin(\omega /2) }{ \omega} \tag{3}$$
and
$$ G(\omega) = \left( \frac{ \sin(\omega/4) }{ \omega} \right)^2 \tag{4}$$
Now, the input spectrum contains all frequencies except the zeros of the sine function which occur at $\omega = 2\pi k$, whereas the output spectrum contains all frequencies except those $\omega = 4\pi k$. Therefore, RBJ is correct in his claim that this particular input-output pair violates the LTI principle of output cannot contain frequencies that do not exist at the input.
However, in my opinion, redefining the output triangular-pulse duration as
$$ \hat{g}(t) = \begin{cases} { 2t+1 ~~~,~~~ -1 < t < 0 \\ 1-2t ~~~,~~~ ~~~~~0 < t < 1 }\end{cases} \tag{5}
$$
eliminates the associated problem, as the output spectrum becomes:
$$ \hat{G}(\omega) = \left( \frac{ \sin(\omega/2) }{ \omega} \right)^2 \tag{6}$$
which includes all frequencies except at those $\omega = 2\pi k$. Thus the missing frequencies at the input are also missing at the output.
Hence, I would like to assume that there was a typo at the duration of the triangular-pulse output $g(t)$, and that otherwise the analysis based on LTI remains valid... Thanks to RBJ for the reminder... EDIT-1-ENDS
EDIT-2 :
However, now happens another problem. If the system is to remain LTI (so that subsequent analysis remains valid), then $\hat{g}(t)$ (the redefined $g(t)$) should be used in the remaning progress. But then the observed waveforms in the figures will be different from what's being drawn. Most notably, the step-response obtained in part-b, becomes a flat output (a unit-step with finite length triangular tails), and finally the impulse response found as step-response's derivative will be a single rectangular-pulse at the origin (and one more at positive infinity) and not a periodic waveform.
So I think the problem should be completely redefined to make any sense...
EDIT-2-ENDS
Your system is stated to be causal and LTI. And an example response for a single rectangular-pulse input $s(t)$ is provided to be a single triangular-pulse output $g(t)$:
$$g(t) = \mathcal{T}\{s(t)\} = s(t) \star h(t) $$
where $h(t)$ is the impulse response of the LTI system, and $\star$ denotes convolution.
The followings are easily deduced from the properties of LTI systems:
a-) Observe that $s_1(t)$ is given by a sum of two shifted $s(t)$ pulses side by side as:
$$ s_1(t) = s(t-0.5) + s(t-1.5)$$
then from the LTI property, conclude the output as:
$$g_1(t) = \mathcal{T}\{s_1(t)\} = \mathcal{T}\{s(t-0.5) + s(t-1.5)\} = g(t-0.5) + g(t-1.5) $$
the output for $s_1(t)$ consists of two triangular pulses repeated (and the origin is shifted).
b-) Now, instead of two pulses, consider an infinity of such rectangular pulses shifted to locate side by side. Then the output will also be infinity of triangular pulses located side by side, as given by the figures in the case b:
$$
\epsilon(t) = \sum s(t -0.5-n) \implies h_{\epsilon}(t) = \mathcal{T}\{ \epsilon(t) \}$$
$$
\begin{align}
h_{\epsilon}(t) &= \mathcal{T}\{ \sum s(t -0.5-n) \} \\
&= \sum \mathcal{T}\{s(t -0.5-n) \} \\
&= \sum g(t-0.5-n)
\end{align}
$$
incidentally, the signal $\epsilon(t)$ is a unit-step function, and the corresponding response $h_{\epsilon}(t)$ is known as step-response of an LTI system.
c-) then the impulse response of the LTI system can be found based on the relation of step-reponse to impulse-response. Namely, if $h_{\epsilon}(t)$ is the step-response of an LTI system, then the impulse response $h(t)$ is :
$$ h(t) = \frac{ d h_{\epsilon}(t)} {dt} $$
you can see from the figure in c-) that the impulse response $h(t)$ is indeed the differentiation of the step-response; i.e, the triangular pulses are converted into bipolar rectangular pulses, and they are repeated indefinetely at the rate of triangular pulses in the step-response. Hence, the figure in case-c.
or it could as well be directly verified as the differentitation property of convolution. Actually this property is a consequence of the differentiation property of the convolution operator...