Improving spectrogram resolution in Python?

Here are the matplotlib.specgram parameters

matplotlib.mlab.specgram(x,
                         NFFT=256,
                         Fs=2,
                         detrend=<function detrend_none at 0x1dd6410>,
                         window=<function window_hanning at 0x1e0b1b8>,
                         noverlap=128,
                         pad_to=None,
                         sides='default',
                         scale_by_freq=None)

The parameters provided in the question description need to be converted to comparable mpl.specgram
parameters.  The following is an example of the mapping:

View range: 0-8000Hz            Fs=16000
Window length: 0.005s           NFFT = int(Fs*0.005) = 80
                                noverlap = int(Fs*0.0025) = 40
Dynamic range: 70dB             n/a
Time steps: 1000                n/a
Frequency steps: 250
Window shape: Gaussian          default window is hanning change to gaussian

If you use 8ms you will get a power of 2 FFT (128).  The following is the description of the Praat settings from their website

  View range (Hz) :
  the range of frequencies to display. The standard is 0 Hz at the bottom and 5000 Hz at
  the top. If this maximum frequency is higher than the Nyquist frequency of the Sound
  (which is half its sampling frequency), some values in the spectrogram will be zero,
  and the higher frequencies will be drawn in white. You can see this if you record a
  Sound at 44100 Hz and set the view range from 0 Hz to 25000 Hz.

  Window length :
  the duration of the analysis window. If this is 0.005 seconds (the standard),
  Praat uses for each frame the part of the sound that lies between 0.0025 seconds
  before and 0.0025 seconds after the centre of that frame (for Gaussian windows,
  Praat actually uses a bit more than that). The window length determines the
  bandwidth of the spectral analysis, i.e. the width of the horizontal line in
  the spectrogram of a pure sine wave (see below). For a Gaussian window, the -3 dB
  bandwidth is 2*sqrt(6*ln(2))/(π*Window length), or 1.2982804 / Window length.
  To get a broad-band' spectrogram (bandwidth 260 Hz), keep the standard window
  length of 5 ms; to get anarrow-band' spectrogram (bandwidth 43 Hz), set it to
  30 ms (0.03 seconds). The other window shapes give slightly different values.

  Dynamic range (dB) :
  All values that are more than Dynamic range dB below the maximum (perhaps after
  dynamic compression, see Advanced spectrogram settings...) will be drawn in white.
  Values in-between have appropriate shades of grey. Thus, if the highest peak in
  the spectrogram has a height of 30 dB/Hz, and the dynamic range is 50 dB (which
  is the standard value), then values below -20 dB/Hz will be drawn in white, and
  values between -20 dB/Hz and 30 dB/Hz will be drawn in various shades of grey.

Link to Praat settings

The OP's question might be concerning the contrast difference between the Praat
specgram and the mpl (matplotlib) specgram.  Praat has a Dynamic Range setting
which affects the contrast.  The mpl function does not have a similar setting/parameter.  The mpl.specgram does return the 2D array of power levels (the spectrogram) the
dynamic range could be applied to the return array and re-plotted.

The following is a code snippet to create the plots below.  The example is ~1m15s speech with a chirp from 20Hz-8000Hz.

import numpy
import pylab
import wave
import array
pylab.close('all')
w1 = wave.open('example_no_noise.wav')
w2 = wave.open('example_noise.wav')
# hmmm, probably a better way to do this, scipy.io function?
x1 = numpy.array(array.array('h', w1.readframes(w1.getnframes())))
x2 = numpy.array(array.array('h', w2.readframes(w2.getnframes())))
x1 = x1 / (2.**(16-1))  # normalize
x2 = x2 / (2.**(16-1))  # normalize
Fs = 16000.
NFFT = int(Fs*0.005)  # 5ms window
noverlap = int(Fs*0.0025)
pylab.figure(1)
pylab.specgram(x1, NFFT=NFFT, Fs=Fs, noverlap=noverlap,
               cmap=pylab.get_cmap('Greys'))
pylab.title('Full 1m15s example min noise')
pylab.figure(2)
pylab.specgram(x2, NFFT=NFFT, Fs=Fs, noverlap=noverlap,
               cmap=pylab.get_cmap('Greys'))
pylab.title('Full 1m15s example more noise')
pylab.figure(3); n=2100*176;
pylab.specgram(x2[n:n+256*256], NFFT=NFFT, Fs=Fs, noverlap=noverlap,
               cmap=pylab.get_cmap('Greys'))
pylab.title('Full ~4s example min noise')
pylab.figure(4); pylab.plot(x1[n:n+256*256])


Conditions for symmetric and unimodal windows in both time and frequency domains

A sufficient condition is, that the distribution is spanned by zero-mean Gauss distributions with non-negative amplitude weights. The spanned space contains only symmetric and unimodal distributions and is invariant under a Fourier transform.


Maximum Ratio Combining vs. Selection Combining

Very quick explanation: SC just throws away $N-1$ of $N$ observations, whereas MRC combines all $N$.

MRC does to a signal in spatial domain what a matched filter does in frequency domain: it maximizes the inner product of signal vector and weights. The exact same derivation as for the matched filter's optimal performance under AWGN apply to the MRC under the same assumption.

You can also do the same consideration on a basic stochastic level.

Assume the combiner has correct knowledge of each observation's SNR.

Then, SC picks the best one. Derivation is simple: if channel realizations are Rayleigh distributed, then the best of these channels has $\sum_{n=1}^N \frac1N$ increased SNR in expectation, $N$ being the number of concurrent observations.

For MRC, you literally get the weighted average of the received signals, so your overall SNR literally is $\sum_{n=1}^N \frac {S_k}{N_k}$, which yields a diversity gain of $N$ (the maximum possible).

Regarding your "I've read papers": the 1959 Brennan paper is the only one you need to read for this :) But, honestly, this is knowledge that's 60 years old, so you don't have to read the original paper, and should probably rather read a good modern textbook or lecture slides explaining combining methods, as these tend to didactically process the paper's content to make it more accessible. A really quick google for "diversity gain MRC" instantly yielded a link to these lecture notes, and you could have found them, too! They're easy to read, I find after quickly scanning them.


Calculating an impulse response in a stable environment

The easiest way to calculate the impulse response is to divide the spectrum of the measured signal by the spectrum of the excitation and do an inverse FFT of the quotient. That's mathematically equivalent to the cross-correlation method but typically easier. It requires that the excitation signal has a reasonable amount of energy at all frequencies so you don't get any "divide by 0" problems.

The length of the impulse response is mostly determined by the reverberation time of the room. For most of home environments that's well below 0.5 seconds so 1s excitation should be plenty. For convenience I would chose a power of two, which makes the FFTs easier. You can still truncate after the IR has been calculated.

Most acoustic systems are actually strict LTI (if no one moves about), however the signal to noise ratio (SNR) tends to be poor especially at low frequencies. A good way to deal with this is to do "coherent averaging". Create a periodic excitation (just keep repeating the sweep), throw away the first two frame acquisitions (long story) and then keep averaging the acquired frames. 10 averages will give you an SNR gain of 10 dB.

Getting the excitation level right is critical: if it's too low the SNR is bad and if it's too high, something will be overdriven and become non-linear. The speakers are probably the first candidate for distortion.


How to evaluate fixed-point implementation of LMS filter is correct?

I'm not an expert on the LMS algorithm. Perhaps you should add a link to the algorithm description so we can help you. However, I have adapted a lot of algorithms to fixed-point implementation, so I think I can help you.
I'm gonna give you a few pointers :
1 - If your algorithm is not recursive, quantization errors can affect the accuracy of your results but the fixed-point implementation will always be stable.
2 - If your algorithm is recursive, it is possible that quantization will radically change the behaviour of your algorithm, it could become unstable, you could be stuck with a limit cycle, etc.
3 - Create a "perfect" model. In your case, the floating-point implementation could be your perfect model. Then, convert your model to fixed-point and compare the results of your fixed-point model. You can toy around with the number of bits for your coefficients, rounding versus truncation, etc.
Regarding your algorithm, you have 6 coefficients, 3 of them are significant (coefficients 3, 4 and 6) while 3 of them seem insignificant (1,2,5). Is this ok?


There are lots of different ways. You could simply put all the image pixel values into a vector with dimension = width x height. :) What you really want to know is how can we describe an image using a lower-dimensional representation and still caption all the important information.

The same approaches that describe an image patch, e.g. SIFT, apply. Just think of the entire image as the patch. However, the method very much depends on the information in the image and the purpose of the vector. For example, you could use some histograms of the gradients, subsample the frequency spectrum, take a certain number of wavelet coefficients, perform edge detection, use moments of the image, divide it into super-pixels and so on. You might need the representation to have different invariances such brightness, rotation, contrast, phase or scale. Basically any method that gives you a condensed version of the information in an image can be put into a vector.

Research into the right method to use for your particular application is probably in the literature somewhere. As for saving a lot of space, this is what compression is all about. e.g. JPEG discards the smallest coefficients of the DCT of blocks in the image to reduce space.


I just answered this in Stack Overflow.  Here's the answer again:

The PNG pixels are in RGBA order, not ARGB, so you'd need to write

alpha | red<<24 | green << 16 | blue <<8

But you are writing color type 2, so your pixels should be 3 bytes each instead of four; you can't really encode them in 4-byte integers as you've done. So either change the color type to 6 or switch to writing each sample as an individual byte.

Also you need a filter byte at the beginning of each row. Zero will work for you. So in your example 4x4 RGB image you will need to write 13 bytes per row

0 R G B R G B R G B R G B
0 R G B R G B R G B R G B
0 R G B R G B R G B R G B
0 R G B R G B R G B R G B

and zlib-compress that. It can all go into one IDAT chunk. If you need to write smaller IDAT chunks, you have to zlib-compress the image first, then split the zlib output into pieces that you put in consecutive IDAT chunks.


A minimum-phase system is always causal and stable by definition, so in the case of discrete-time systems with rational transfer functions, all poles and zeros are inside the unit circle of the complex $z$-plane. This is why a minimum-phase system can be inverted by a causal and stable system. Because of this definition, it does not make sense to compare non-causal systems and their phase lag or group delay to the corresponding properties of minimum-phase systems. Minimum-phase systems have the smallest possible phase lag of all causal systems with the same magnitude response.

The group delay is a function of frequency. For symmetric/anti-symmetric impulse responses, the group delay is constant. This is probably what you meant in your question. Such linear-phase (constant group delay) systems can never be minimum-phase (unless their impulse response is just a unit impulse at $n=0$). Minimum-phase systems have a frequency dependent group delay, but at any frequency, their group delay is smaller than the group delay of any other causal system with the same magnitude response.


What Are the Limitations of a Canny Edge Detector?

From my experience, the following points are limitations:

The result is binary. You sometimes need a measure of 'how much' the edge qualifies as an edge (e.g. intensity image coming from a Sobel amplitude edge detector)
The amount of parameters leads to infinitely tweaking for getting just that little better result.
You still need to connect the resulting edges to extract the complete edges that seem so obvious for the human eye+mind.
Also due to the gaussian smoothing: the location of the edges might be off, depending on the size of the gaussian kernel.
The method has problems with corners and junctions:

The gaussian smoothing blurs them out, making them harder to detect (same goes for the edges themself)
The corner pixels look in the wrong directions for their neighbors, leaving open ended edges, and missing junctions

This last problem is addressed by the SUSAN method, which connects edges better and also results in nice junctions, as shown by these example figures as given in the linked paper:

Test input image:

Results SUSAN:

Results Canny:

You can clearly see SUSAN finds the corners and junctions instead of Canny.


How to average MFCCs from different speakers saying the same word?

What you can try is the DTW (Dynamic Time Warping) distance (https://en.wikipedia.org/wiki/Dynamic_time_warping).

for each instance of the spoken word "west", you compute MFCC's per frame and then calculate the DTW distance between each of the pre-recorded instances of word. this gives the rough estimate of the valid DTW distance. now when a new instance of the word comes we again calculate DTW with one of the earlier instance and if the DTW distance lies roughly in the range calculated before, it is the same word.


Inverse DFT to a real signal when only odd harmonics are present

So, I've read the paper I cited more attentively, and it turned out that the author, while not having described the inverse algorithm explicitly, left enough clues in the derivation of the forward algorithm, after all. In case anyone needs it in the future, here it is.

Input: the nonzero (odd) harmonics $X_k$, $k=0,\dots,N/4-1$ (assuming here for simplicity that $N$ is divisible by 4).

Build the array $Y_k$, $k=0,\dots,N/4$: $Y_0 = -2\mathrm{Im}X_0$, $Y_{N/4} = 2\mathrm{Im}X_{N/4-1}$, for the rest $Y_k = i(X_{k+1}-X_k)$.
Calculate $y = \mathrm{IRFFT}(Y)$ (the standard inverse FFT to real signal, e.g. numpy.fft.irfft in Python or dft_c2r from FFTW). The resulting vector $y$ is real and has length $N/2$.
$x_0 = \frac{2}{N}\sum_{k=0}^{N/4-1} \mathrm{Re}X_k$, $x_n = y_n / (4 \sin(2\pi n))$, $n=1,\dots,N/2-1$. This gives the first half of the full vector $x$.

Of course, this method is not ideal, since the division by a sine can introduce instabilities. There may exist a method to deal with that (there is one for the FFT of an (anti)symmetric signal, which requires a similar post-processing).


The answer to this question could be tricky.
For a simple answer you may take the bandwidth of a bandpass signal as $\omega_2 - \omega_1$ (highest minus the lowest frequencies), and don't care about what happens in between the two. Hence your bandwidth allocation is not fully efficient, and a sampling based on this bandwidth will include redundant information of those empty bands.
Thus, for a more rigorous answer about finding the absolute minimum number of samples to represent a bandpass signal whose spectrum may include empty bands (gaps) between the minimum and maximum fequencies, you would take only the nonzero intervals into count.
Yet when the signal is to be transmitted over a physical channel without frequency allocation control, then you would again consider the minimum and maximum frequencies as the bandwidth, unless cross channel spectrum ovarlapping is allowed.
So the answer depends on your intentions.


How does the number of frequency-domain sampling points influence the outcome of an inverse FFT?

Regarding example 1:
first of all, either the fft or the ifft needs to be normalised by the number of sampling points as you have done (actually you can normalise each by a factor of the square root of the number of points, it is just a meter of definition). However, in your case the ifft is half the length of the fft you performed. Hence your signal is divided by 2. If you would perform the normalisation when calling the ifft (using half the original number of samples) you should get the correct height.
Now, regrading the second part of the signal, it is not truncated. What you see is some kind of aliasing. Since you downsampled  without filtering, you get aliasing. This is the same concept as frequency aliasing (DFT and IDFT are basically the same). You should try to flip the input signal in the time domain and witness the results.

Next, regarding example 2:
Your problem is with the number of sampling points. Let me start by mentioning that since your signal is real, the frequency domain is symmetric.
In your first example you had even number of points. when you downsampled the frequency you still got a symmetric signal. In the second example you have odd number of samples. Therefore your decimated frequencies are not symmetric and the result is a complex signal. If you try it with one less (or more) sample it should give a similar result to example 1. On the other hand you can perform the frequency decimation while keeping the signal symmetric.


Convolve a minimum phase wavelet with delta peak

If I understand the question properly. Consider the case when the delta is at time zero. In this case the convolution does not do anything to the other signal i.e. it starts at time y and the maximum peak will be at time z.

Now, if we delay the delta to time x, then the output of the convolution will also be delayed by x. So the signal will start at time x+y, and the maximum peak will be at time x+z.

If I haven't understood correctly, you may need to include some figures to more clearly explain the situation.


Consider a discrete-time input signal of the form:
$$ x[n] = \cos(\omega_0 n) ~~~,~~~-\infty < n < \infty, ~~~~~ n\in \mathcal{Z}$$
where the radian frequency $\omega_0$ is set between 0 and $\pi$ radians per sample.
Now, consider two basic discrete-time (digital) filters which are defined through addition and subtraction of their input $x[n]$ for producing their outputs $y_1[n]$ and $y_2[n]$ as:
$$y_1[n] = (x[n]+x[n-1])/2 ~~~,~~~ \text{sum filter}$$ and
$$y_2[n] = (x[n]-x[n-1])/2 ~~~,~~~ \text{difference filter}$$
Lets make a qualitative analysis of these two filters by setting their input frequency $\omega_0$ to low (close to $0$), and high (close to $\pi$) values, and then observing the corresponding outputs respectively;
First, assume that $\omega_0$ is set to low frequencies: Then the consecutive input samples $x[n]$ and $x[n-1]$ will have highly similar values, as a low frequency sine wave will not change much from one sample to the other. When this is the case, their sum will add up, whereas their difference will cancel.
Therefore, $y_1[n]$ will be approximately equal to the input's value $x[n]$, while the output $y_2[n]$ will be close to zero, due to cancelling of $x[n]$ by $x[n-1]$ when subtracted. Then we conclude that the sum filter passes the low frequencies while the difference filter attenuates (blocks) them.
For the second part of the analysis, set $\omega_0$ to high frequencies: then the values of the input samples $x[n]$ and $x[n-1]$ will be of opposite polarities, as the cosine will be rapidly changing from sample to sample. Then, their sum  will cancel, whereas their difference will add up.
Therefore, $y_1[n]$ will be approximately zero, while $y_2[n]$ will be following the input. Then we see that the sum filter blocks high frequencies while the difference filter passes them.
Combining these two analyses we conclude that the sum filter is a lowpass filter, which passes low frequencies and attenuates high frequencies, while the difference filter is a highpass filter, which attenuates low frequencies and passes high frequencies.
Following this basic basic setting, more complex filters are realized through using more samples at farther delays, and are weighted properly. Then the passband and stopband cutoff frequencies, the transition bandwidth and ripples are all determined by those delays and weights applied into the summing (or differencing) samples and the number of samples (length of filter) used in the sums (or differences).
Those weights are then called as the filter coefficients (or its impulse response $h[n]$ for the FIR filter) that characterize the filter.


Matlab does not accept lesser than $1$ indices. Other languages often accept a $0$-based index. More about Shifting Sequences Using MATLAB is available, including function circshift.m.


Question about harmonic calculations and other functionals for image analysis

Perhaps looking at this paper 1 would be appropriate for a literature survey.

The paper describes in detail the rules of constructing functionals for a Trace transform. The functionals described in your paper are merely examples of the functions which follow the rules presented in paper 1.


Frequency dependent phase difference between signals arriving at two microphones

For a constant time difference, the unwrapped phase difference measured at the centers of the two windows is linearly related to frequency. (To measure phase at the center of an FFT window, use an fftshift prior to the FFT, or flip the phase of the odd numbered FFT result bins.)  So you can use all of the FFT results, not just one bin.   The low frequency bins may help in determining the potentially ambiguous phase unwrapping of the higher frequency bins, if the signal is sufficiently wideband.


Could someone please explain how to get frequency response of a given Bode Plot?

If your question is how to use this graph to find the transfer function the answer is simple : you can't.

First of all, this graph give only half of the information: you have a gain plot but not the phase plot, and 2 systems with different TF could have the same gain profile.

Second of all, this is obviously an asymptotic graph, as it presents a discontinuity around the cutoff frequency.

However, you could extract some information from this plot. You can see that it cuts the high frequencies with a slope of -20dB per decade, which means that there is a great chance that it is a 1st order low pass filter. Additionally, you can also read that  the cutoff pulsation is about 0.1 rad/s.

Think of it like a normal function, would you be able to tell what a function is if you any saw an asymptotic graph of one half of it? All you can do is make some educated guesses about the form of the TF.


What is Causing a Noise Signal in USB Audio Recording?

It's NOT USB problem. If you there was anything wrong with the digital connection, you would get drop outs or nothing at all. The error detection/correction on the USB bus is binary: either you got all the bits or you didn't.

This is clearly analog noise, probably from switching power supply. It happens most likely in the analog front end of the mixer. Once it's digitized, it can't pick up this type of noise anymore.

You can try a different USB recording device or throw a USB hub in the path to isolate the iMAC power from the Xenyx power see if it helps. It's also entirely possible that the Xenyx will do this no matter what: Behringer has a reputation for great value, but also for sloppy design and poor quality control. The Xenyx is a blatant ripoff of the Mackie Onyx (as the name implies).


How to quantify area of interest in an image?

Given that the size of the barnacles is not expected to be considerably larger than the area of the hull being imaged at a certain shot, you can simply deduce the actual barnacle coverage area by counting how many pixels belong to the hull and then multiplying by a "pixel size" that you can derive by including an item of known length in the image.

The key assumption behind this (simple) technique is that the barnacles and the hull exist on the same plane (As if the barnacle layer was a "sticker" on the hull).

Obviously, this is not the case because the barnacles are three dimensional and "stick out" of the hull. For the same reason, the angle at which a certain part of the hull is imaged might also introduce an additional error factor.

You can control for all this by imaging relatively large areas of the hull (i.e. the camera is at a relatively long distance from the hull), taking care for the camera to be pointing vertically to as much of the hull as it is possible (this of course conflicts with the "take a picture at a long distance" and some experimentation will be required). Finally, you will also need to have good control of the lighting of the shot and for this reason, you might have to use a flash bulb (even in day-light) to flood the scene uniformly with light and avoid shadows.

Having taken the image, you can then convert its colour space from Red, Green, Blue (RGB) to Hue, Saturation, Brightness (HSB) and examine the Hue component's histogram for the uniform colour of the hull (here, Red and because of this, relatively easy to pick out). Again, this assumes that the hull is painted using just one colour.

Here is the process in Gimp:

The three small B&W images at the botoom of the screenshot are the HSB components respectively. The histogram at the top right is derived from the H component (the first B&W image at the bottom left-hand side of the screenshot).

As you can see, there is a clear peak towards the beginning of the H spectrum corresponding to the red colour (Please see the H band in this image to understand how the horizontal axis of the histogram now "corresponds" to hue value).

If there was an object of known length on the hull (an iron waterline marking relief for exaple or a pen or a coin) we could derive the pixel size and through this, the actual surface of the whole image. By then counting how many pixels "belong" to the peak we could deduce the actual surface occupied by the barnacles.

For example:

Image Dimensions: $450 \cdot 229 = 103050$ pixels

Pixels corresponding to red: $8494$

(Please note:You probably cannot see this because of the size of the image but once you have your histogram in GIMP (Colours->Info->Histogram) you can mark an area on the histogram and read how many pixels it describes. This is how this number was derived.).

Barnacle coverage area (%): $\frac{103050-8494}{103050} * 100 \approx 91.7\% $

All you need now is an object of known size to convert this to $cm^2$

In terms of automating this process, it can be done using a number of techniques and platforms (e.g. Python, Octave, etc) but for a more accurate answer reagrding those, a little bit more detail would be required regarding the application.

Hope this helps.


Which is the Matlab equivalent to IDL's function bytscl?

Here is my answer, thanks for the comments, this version may be wrong but anyway the image data looks better.

function [ Y ] = bytscl( X )

m = min( X( : ) );
M = max( X( : ) );
Top = 255;

    if isfloat( X )
    Fc = ( Top + 0.9999 )*( X - m )/( M - m );
    Y = cast( Fc, 'int8' );
    else
    Ic = ( Top + 1 )*( X - m - 1 )/( M - m );
    Y = cast( Ic, 'int8' );
    end

end


$E\{|D_k|^2\}$ is the expected value of the squared magnitude of the noise in the $k^{th}$ frequency bin. Since they assume that the noise is zero mean, this equals the noise variance. Remember that the variance of a random variable $X$ is defined as $\sigma_X^2=E\{|X-\mu_X|^2\}$, where $\mu_X=E\{X\}$ is the mean of $X$. So if $\mu_X=0$ (as assumed here for the noise) then $\sigma^2_X=E\{|X|^2\}$. Note that $E\{|D_k|^2\}$ is also the noise power in the $k^{th}$ bin, so it can't be zero unless there is no noise at all in that bin.


Get I and Q components from a signal (multple frequency signal)?

If you need to do it in time domain you need an ideal Hilbert transform filter and multiply it sine function, and pass it through low pass filter (lpf)to get Quadrature (Q) phase. but don't expect accurate results, or maybe you can use PLL but you have to use it brilliantly like below.


Is This Interpretation of Auxiliary Independent Vector Analysis (AuxIVA) Correct?

I think your understanding is right, and it has already been implemented in many open sourced codes, like https://github.com/LCAV/pyroomacoustics and https://github.com/n-narisawa/bss

However, the performance of AUXIVA is not as good as original IVA according to my experiements. Do you have some new conclusions about the performance and the real-time implementation of AUXIVA? Thanks.


intuitive interpretation Fourier Transform of two different Rectangular pulses

Note that the peak height of the Fourier transform of a wider rectangular function is proportionally higher.  This is necessary to preserve total energy, as per Parseval's theorem.  However, the energy required to create the sharp edges of a rectangle in the time domain is still out there on the tails of the Sinc function in the frequency domain, roughly the same (but more wiggly).

But when you scale the vertical axis of 2 plots so that they have the same peak height, the envelope of the FT of a wider rect function looks narrower in the plots of that transform.  But that's only in relation to the peak with a scaling applied.


Linear prediction (LPC) of Sine wave samples around maximas

Disclaimer: This is explanation is based on observation of my MatLab plots and my note be 100% textbook correct!

Well after much looking around and experimenting I read something about the 'yule-walker' method for estimating the coefficients assuming the signal to be 'zero' outside of the observation point. This is not the case for my co-/sine-wave, so the estimated coefficients and the resulting predicted samples $\hat{x}(n)$ will trend towards zero. The further away from the $x$-Axis the sharper the slope towards 'zero' when using yule-walker.

For my problem an approach using Burg's method (coeffs = arburg(x_kn, P);) seems to yield the result I'm looking for, since burg does 'not' assume the Signal trending towards 'zero'?? At least this is how I can explain the result of my two plots. The 'star' marks the point where the prediction of future values start.

Yule-Walker:

Burg-Method:


I want to find the first period $T_0$ in a signal $x[n]$ and save it to a vector?

The time period should be inverse of the signal's frequency

  T_0 = 1/freq

not sampling frequency


I need to insert an underwater acoustic channel model to an existing code that compared between OFDM and GFDM

I modeled the the acoustic channel noise in matlab:

    % On the Relationship Between Capacity and Distance in an
    % Underwater Acoustic Communication Channel

    Turbulence=17-(30*log10(f));
    Noise_Turbulence=  power(10,(Turbulence*0.1));

    Shipping=40+(20*(s-0.5))+(26*log10(f))-(60*log10((f)+0.03));
    Noise_Shipping= power(10,(Shipping*0.1));

    Wind=50+(7.5*(power(w,0.5)))+(20*log10(f))-(40*log10((f)+0.4));
    Noise_Wind=power(10,(Wind*0.1));

    Thermal=-15+(20*log10(f));
    Noise_Thermal= power(10,(Thermal*0.1));

    Noise=Noise_Turbulence+Noise_Shipping+Noise_Wind+Noise_Thermal;
Ambient_Noise=10*log10(Noise); %the values will be converted into dB


Why is there only one integration in the solution if there is two integral in the formula?

You have a single random variable $\Theta$ in this example, so taking the expectation with respect to that random variable results in a single integral. The formula for $R_X(\tau)$ in your question is rather messy and full of mistakes; it should be as follows:

$$\begin{align}R_X(\tau)&=E[X(t+\tau)X(t)]\\&=E[a\sin(\omega_0t+\omega_0\tau+\Theta)\cdot a\sin(\omega_0t+\Theta)]\\&=\frac{a^2}{2}E[\cos(\omega_0\tau)-\cos(2\omega_0t+\omega_0\tau+2\Theta)]\\&=\frac{a^2}{2}\frac{1}{2\pi}\int_{0}^{2\pi}[\cos(\omega_0\tau)-\cos(2\omega_0t+\omega_0\tau+2\Theta)]d\Theta\\&=\frac{a^2}{2}\cos(\omega_0\tau)\end{align}$$


Attenuate the low frequencies in a speech wave

As you've already figured out, you need to choose an odd filter length. I'll explain why. You're designing a linear phase FIR filter, and such filters satisfy certain symmetry conditions. See this answer explaining the four types of linear phase FIR filters (even/odd length, even/odd symmetry).

For frequency selective filters, such as the one you're trying to design, you need filters with even symmetry (type I or type II). An even length (type II) filter necessarily has a zero at Nyquist, so you can't use it to design a high pass filter. A type I filter (even symmetry, odd filter length) doesn't have this restriction, so in order to design a linear phase high pass filter you need to choose an odd filter length, i.e., a type I linear phase FIR filter.


You are correct that if your signal is bandlimited to <5 Hz, then you can perfectly represent it with a 10Hz sampling rate.  This is the well-known sampling theorem

But ... there may be practical considerations for why one would not be able and/or inclined to use critically sampled data.

One reason is the difficulty of making a signal critically sampled. Any operation you perform to change the rate of the signal is going to have some filter with a non-zero transition bandwidth.  In your example, this limits the unaliased frequency content to 5-ftrans  This transition bandwidth can be made very narrow with long impulse response filters but this has costs both in terms of processing and in transients (ringing) at signal start and end.

Another reason is the efficacy of algorithms that work on the resulting signal. If you need to work with a blackbox component that can only choose the nearest sample, then you'll be better off feeding it oversampled data.

Most (all?) non-linear operations will behave differently with critically sampled vs oversampled data.  One example is squaring a signal, a well known method of BPSK carrier recovery. Without a 2x oversampled condition, the multiplication of the time domain signal with itself causes wraparound garbage aliasing when the frequency domain convolves with itself.


What are some ways to detect the fundamental frequency in a signal using without using FFT?

Count the time interval between zero-crossings (if the waveform/SNR allows...).
Autocorrelation, average magnitude difference function and their variants (YIN...).
Estimating the parameters of a damped sinusoids model (MUSIC/ESPRIT).
Several instances of Goertzel algorithm, or more generally other filterbanks (band-pass, comb) - in case there is only a small number of candidate frequencies that you want to discriminate.


Using DFT to determine control system bandwidth from user input

Your basic approach seems correct! However, step 2 is not a good idea since taking the square of the signal is a non-linear operation. Taking the square of the FFT to determine the signal power is correct though.

For your profs. question it depends on what you have measured and at which sample rate. I assume it is the steering wheel angle and therefore it should be possible to answer it, but there will be some pitfalls as most likely the steering movement is not in a sinusoidal form. However, the basic frequencies in the movement will be related to this.

Given that you redo the estimation without the squaring of the signal, the ten times safety margin should be more than enough.

To further convince you and your prof. I would suggest bandlimiting the measured data to your estimate by using a very sharp FIR filter and compare visually as well as downsampling the data and compare.


What Noise Will Change the Frequencies Components of a Signal?

I'm not sure I understand your question.

But usually the noise model is additive.
Namely, due to the linearity of the DFT, at each bin you will have the summation of the signal contribution and the noise contribution at this specific frequency bin.

Now, let's say we have 2 bins.
At the first the signal is 7 at the second is 6.
The noise is 3 and 5 respectively, now the summation will shift the DFT peak though the signal isn't changed.


Measure of intermittency/continuousness of a signal

Any metric measuring the relative importance of high-frequency content vs low-frequency content will be correlated with the amount of abrupt changes in the signal.

You can try spectral roll-off (the frequency under which x% of the energy is concentrated) or spectral centroid.


The fact that you are plotting the FFT of the whole audio clip makes me think that you are looking at the wrong class of solutions. From the waveforms of signal A and B, it is clear that your noise and signal are not stationary. If you want to use a frequency-domain noise reduction technique, this should be done, in your case, on shorter overlapping windows across which you can assume that the properties of your signal and noise are stationary. You cannot avoid using the Short-term Fourier Transform.

In the rest of this post, I will switch notations to:

$X$ original (not noisy) signal
$N$ noise signal
$Y = X + N$ your observation, that is to say $B$. $X$ is what you want to recover, $N$ in your case is known (signal $A$)

Most denoising methods share the same structure: get a STFT representation of $Y$ and $N$ (matrices $Y(p, k)$ and $N(p, k)$, where $p$ is the frame index and $k$ the frequency bin index), process pairs of FFT frames from rows of $Y$ and $N$ to get an estimate $\hat{X}$ of a slice of the denoised signal spectrum, and use it to build a filter/mask applied in the frequency domain to $Y$, the conversion back to the time-domain being performed through overlap-add.

I suggest you to try to implement the Ephraim-Malah algorithm, which is a classic and often cited "baseline" denoising method. It is based on two intuitive tricks:

Estimate the SNR and attenuate the strength of the denoising filter in segments with low SNR (this prevents many artifacts in areas where the signal is close to the noise level).
Use temporal smoothing (if we have a high SNR in a frame, make the filtering/subtraction stronger in the next frame since it is also likely to contain signal...).

I have found a relatively self-contained set of slides on the topic - the author introduces increasingly sophisticated denoising rules, from spectral subtraction to full-blown Ephraim-Malah. The only thing he doesn't tell you to get to a practical implementation is how to get the estimate $\hat{S}_x(p,k)$ - since you don't observe $X(p,k)$, just $Y(p,k)$ and $N(p,k)$. The simplest approach is to subtract the power spectrum of the noise from the power spectrum of the signal, and zero the STFT cells in which the result is negative:

$\hat{S}_x(p,k) = \begin{cases} |Y(p,k)|^2 - |N(p,k)|^2 & when |Y(p,k)|^2 - |N(p,k)|^2 > 0 \\ 0 & otherwise \end{cases}$

You can then try the filtering rules shown on pages 5 and 10 of the slides.


People have been decoding FSK signals for ages with very weak microprocessors.  The usual method is called zero crossing detection.

You basically just measure time (count the number of samples) between the times when the signal crosses the zero volt line.
The number of samples multiplied by the sampling rate tells you the rate as a time period.  One divided by that number is the frequency of your signal.

That works well enough with clean signals.  If your signal has noise mixed in with it, then you use a band pass filter to just pass the frequencies you are interested in.

I'd work up an FSK demodulator using zero crossing that works for clean signals.  Once that works, you can worry about filters and improving the performance.


A common technique for computing the inverse fft is to invert the imaginary part of the input array, perform a forward fft and then invert the imaginary part of the output array. In the case of the Cepstrum, the input to the ifft is real-valued (due to the absolute value function) and also symmetric. In this case there is no imaginary part on the input to invert, and the symmetry implies that the ifft will produce no imaginary components on the output that would need to be inverted.  Hence in this case the fft and the ifft are fully equivalent.


Do the DCT coefficients correspond to particular pixels?

What I'm planning to do is to apply edge detection to obtain the
  locations/indices of the pixels and then get their corresponding DCT
  Coefficient

As others stated in the comment - any pixel in the given block is related to all co-efficient in the block and vice-versa.

However, understand that Edge detection is a process of finding the gradient. Linear convolution essentially does the job. What you can do is to apply equivalent of that convolution in the DCT domain which will give you an idea about the edge orientation in the DCT block domain.

While, this is relatively narrow area - there are papers which allows you to find edges in DCT blocks. See the references below:

A Compressed Domain Scheme for Classifying Block Edge Patterns by Hyun Sung Chang, and Kyeongok Kang IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 14, NO. 2, FEBRUARY 2005 145
Direct feature extraction from compressed images by Bo Shen et. al.
SPIE vol. 2670, Storage & Retrieval for Image and Video Databases IV, 1996


Normally if you sample at Nyquist rate it's "good enough" but how does this interact with the discrete time correlator?

It interacts very well! See the following:

  Does sampling at a higher rate benefit the system?

You mean oversampling.

If you filter your signal to the original signals bandwidth, which is now smaller than the Nyquist rate. You thus get a component that is correlated between samples (which is your signal) and one that isn't (which is your noise). By averaging / downsampling your signal, you increase the variance of your signal component more than that of your noise component – you get an SNR gain.

  Does up-sampling the signal and passing it through the correlator yield the same benefit?

No, because upsampling doesn't "uncorrelate" the noise that you already digitize.

  Is there any benefit to over-sampling and then discrete time matched filtering?

Yes, as described above, the oversampling allows you to increase your SNR. Think of it this way:

Assume that if critically sampled, your pulse length is $N$; thus, your MF does a correlation over $N$ samples.

You oversample your analog signal by a factor of $f>1$. You interpolate your MF by the same factor, so that the pulse and the filter are now $Nf$ in length.

While the power of the samples stays the same, the amplitude of the signal part in the  (digital) correlator output has been scaled by a factor of $f$, which means the power  has been scaled by $f^2$.

For the noise component in the correlator output: White Noise samples are uncorrelated, and we know that thus the variance of their sum (which is what a correlator is, weighted sum) is the sum of their variances. Since variance is power for zero-mean noise, that means we get a noise power increase of $f$.

Thus, oversampling by $f$ gets us an SNR increase of $\frac{f^2}f$ :)


Why correlation is not associative in the context of image processing?

Subtraction is also not associative: (a-b)-c does not equal a-(b-c) in general. The order of the flipping of the sign matters, just like it does for the flipping of the kernel.


Zero order hold interpolation and Nearest-neighbor interpolation?

Zero-order hold and nearest-neighbor interpolation are conceptually the same thing; when extending the original dataset to cover new data points, you are just picking a value in the same area of the signal to use as the interpolant, instead of something more sophisticated.
The main difference I could see between the two is that zero-order hold is usually used in the context of interpolation of a 1-D signal, particularly to describe the action of a digital-to-analog converter. In that model, the zero-order hold is usually rendered as starting at the corresponding sample time $n$ and extending until the next sample instant $n + 1$.
Nearest-neighbor interpolation is usually drawn somewhat differently, where the region of influence of a particular sample is centered on the sample position, so instead of the value being held from sample times $n$ to $n + 1$, that sample value would instead influence the interpolated output from $n - 0.5$ to $n + 0.5$.


In fact you have two signals, and it depends on what you want to achieve, but usually you would just filter both signals (the real and the imaginary part) with the same (real-valued) low pass filter. So you either need two (identical) low pass filters, or you filter both signals sequentially with the same filter.

In the general case, filtering a complex signal $x(t)=x_R(t)+jx_I(t)$ with a complex impulse response $h(t)=h_R(t)+jh_I(t)$ requires four real-valued filtering operations (just like complex multiplication requires four real-valued multiplications):

$$\begin{align}y(t)&=(x\star h)(t)\\&=(x_R\star h_R)(t)-(x_I\star h_I)(t)+j\left\{(x_R\star h_I)(t)+(x_I\star h_R)(t)\right\}\end{align}\tag{1}\\$$

In your case $h_I(t)$ is zero, so you're left with just two filtering operations.

One final remark on complex convolution: actually, one can get away with only 3 real-valued convolutions (just like complex multiplication really needs only 3 multiplications if you're smart):

$$\begin{align}y_1(t)&=(x_R\star h_R)(t)\\
y_2(t)&=(x_I\star h_I)(t)\\
y_3(t)&=((x_R+x_I)\star (h_R+h_I))(t)\\
y(t)&=y_1(t)-y_2(t)+j\left\{y_3(t)-y_1(t)-y_2(t)\right\}\end{align}\tag{2}\\$$


Calculating penalty to SNR for non-matched filter

Let's consider a received signal

$$Y(t)=Ap(t)+N(t)\tag{1}$$

where $A$ is the information symbol (modeled as a random variable), $p(t)$ is the transmit pulse, and $N(t)$ is additive white Gaussian noise (AWGN, modeled as a random process). For the sake of simplicity let's assume that all signals are real-valued (baseband case).

Filtering with an LTI system with impulse response $h(t)$ and sampling at $t=0$ gives

$$Y_0=\int_{-\infty}^{\infty}Y(t-\tau)h(\tau)d\tau\Big |_{t=0}=\int_{-\infty}^{\infty}Y(-\tau)h(\tau)d\tau\tag{2}$$

With $(1)$ you get

$$Y_0=A\int_{-\infty}^{\infty}p(-\tau)h(\tau)d\tau+\int_{-\infty}^{\infty}N(-\tau)h(\tau)d\tau\tag{3}$$

The power of the desired signal (the first term in $(3)$) is

$$P_S=E[A^2]\left(\int_{-\infty}^{\infty}p(-\tau)h(\tau)d\tau\right)^2\tag{4}$$

where $E[\cdot]$ is the expected value, and $E[A^2]$ is the power of the information symbol $A$. From the second term in $(3)$, the noise power is

$$P_N=\frac{N_0}{2}\int_{-\infty}^{\infty}h^2(\tau)d\tau\tag{5}$$

where $N_0$ is the one-sided power spectral density of the noise. From $(4)$ and $(5)$ the SNR is

$$SNR=\frac{P_S}{P_N}=\frac{E[A^2]\displaystyle\left(\int_{-\infty}^{\infty}p(-\tau)h(\tau)d\tau\right)^2}{\displaystyle\frac{N_0}{2}\int_{-\infty}^{\infty}h^2(\tau)d\tau}\tag{6}$$

We know that $(6)$ can be maximized by a matched filter:

$$h_{opt}(t)=Cp(-t)\tag{7}$$

with some constant $C$. From $(6)$ and $(7)$, the optimal (maximum) SNR is given by

$$SNR_{opt}=\frac{E[A^2]}{\frac{N_0}{2}}\tag{8}$$

The actual SNR for an arbitrary receive filter impulse response $h(t)$ can be written in terms of $SNR_{opt}$:

$$SNR=SNR_{opt}\underbrace{\frac{\displaystyle\left(\int_{-\infty}^{\infty}p(-\tau)h(\tau)d\tau\right)^2}{\displaystyle\int_{-\infty}^{\infty}h^2(\tau)d\tau}}_{\le 1}\tag{9}$$

where the last term in $(9)$ (or its inverse) could be called the "SNR penalty".


Can I have some feedback on this research paper?

Comments provided here are in two broad categories: Presentation and Subject matter. The "Presentation" section is the easiest to address. There are some things that could be rephrased in terms of language use but these might be just personal preferences. The "Subject matter" section includes comments in methodology which could take more time to address.

Overall, my perception was that the paper is an excellent start to more complex work but certain key points for this type of work were not clear enough. There is always "more work" that can be done but this particular paper has room for improvement before hitting those limits. I would give it another go before submitting to a journal.

Presentation

"Passive noise control uses bulky sound absorbent materials to mas unwanted sound, whose performance...", The control's or the material's?
"In the existing ANC literature, which has a history of 7 decades, the acoustic...", do we really need to load that phrase with yet another point on the way to saying something about the existing literature? It breaks the flow and I am not sure how much it adds to the point.
"This approximation, used for ease of computation, limits researchers....". Same comment applies here. Just say what you want to say in one phrase with two strong points rather than having one main point and an aside.
"...treatment of the research work is available at [some URL]...". Why not reference this properly and include one more citation? Or add a footnote if you think it would be better to have the URL immediately visible for some reason. As a reader, I am not even going to "parse" the URL unless I want to visit it and I have not yet read the message that this paper is trying to convey. Why send the reader away already?
"...using single or multiple microphones and passes [IT] on to the...","...THE TOOLBOX provided by Habets[20] HAVE been..." I am not going to flag each and every one of these tiny little things, please proof read the document. A typo can be seen a bit like "Oh come on man! I have worked my bottom off in this paper and you are picking on a typo!?!?!". But the thing is, the typo is brakng the flw. You are nt the redr of the papr, someone else is and they have no idea what you are trying to tell them. They "form" the idea by reading these lines, word-by-word. Typos get in the way. I stop, re-read the phrase, start again...It's more effort.
"ANC PRELIMINARIES" goes on for 1.5 paragraphs BEFORE revealing that what is described is actually depicted in Fig 2. For 1.5 paragraphs, I have been building this sketch up in my head reading the description. More effort on my part (the reader...not me personally).
"...Room Acoustics Models(RAM) are often used to recreate and reproduce room acoustics without any cross-model mismatch..." What is "cross-model mismatch"?
"...beamforming, acoustic auralization, speech.." Auralization is already acoustic.
"A RIR model channel [Fig 3]...". The figure's caption could have RIR captialised "Room Impulse Response...".
Are you sure you need II and III as separate sections? I don't know the audience of the paper, but maybe these two sections are at a kind of...very introductory level (?).
You start IV with "Figure 4 shows...", that's great. Why not do the same with the rest of the figure references?
Caption of Table 1 reads "Simulation settings used for an arbitrary case". We will come back to the "arbitrary case" later on in section 2 but as far as presentation of the work is concerned "ARBITRARY CASE" reads oddly. Perhaps this is the "TEST CASE" or "SIMULATION PARAMETERS"?
Add a grid to all of your graphs. It helps with reading the graph a lot.
"The source signal considered is sinusoidal with a wideband component", "We used the same method to generate the source signal as used here[22]...". Everything that is key to the results of this paper must be sufficiently mentioned in THIS paper. I (the reader...) am not going to go to your github repository to check your code. I may have a system of my own and I am just looking to improve or test it. I would really appreciate a simple formula there showing to me how was the "sinusoidal with a wideband component" created.

Subject Matter

Section IV:

"The room size was chosen arbitrarily to a size of...". Why arbitrarily? What is wrong with "The average auditorium size is...therefore we focused on rooms of that size", or, "In our research we chose to focus on rooms of size BLAH because that's the standard listening room size". Maybe you can refer to some sort of authority that is sending standards for room sizes. Maybe the equipment you develop is ideal for rooms of BLAH size? In any case, for me "...chosen arbitrarily..." is not perceived positively. The point is why did you choose that room size?
"The z level of [THE] microphone, noise source and microphone were fixed to a height of 1.53..." WHY?. First of all, 1.53 (what? Probably meters) and then WHY?. Maybe that's the level an average person's ear ends up when they are sitting? Otherwise...WHY 1.53 meters? Is that the average size person? Is there any other reason?
"...for ease of computation....". No.
"The x coordinate[S] [...] were chosen randomly...". Why randomly? What is the reason? I don't know if it is clear what I am saying here. I am not trying to put you on the spot. What I am reading is "We took some tools, we threw some dice, we got these numbers"...Which, I don't think is what you are trying to say here. Pin these parameters down because of a reason. Then the paper gets some perspective with reality.
"The antinoise source positions were evenly spaced at a distance of 11cm". WHY?.
"Table 1 shows an arbitrary case..."... Shows one of the test cases considered, see previous comment... There is too much "arbitrary", "random".
*"The sampling frequency (Fs) is fixed to 2kHz"....WHY?. For example, you could say something like, "Given the current performance of DSP / speakers, realistic ANC is only possible for frequencies between BLAH and BLAH Hz. We targeted this band and therefore worked with Fs=2kHz.".
"...according to Nyquist theorem...". No.
*"The adaptive filter [ORDER] was fixed [AT] 350 [COEFFICIENTS] at Fs=2kHz which is an important parameter determining ANC performance."... WHY?. Why is 350 at 2kHz important for performance? Are you talking about computational performance or whether the ANC does its job?
"A step size of .." It is not entirely clear that you "sweeped" the locations of the mics and the speakers.

THE POINT ABOUT FIR FILTERS: In the introduction it is stated "In the existing ANC literature, [...], the acoustic channels are mostly modelled using [FIR] filters. This approximation, [...], LIMITS researchers in studying the effects of physical aspects like transducer positioning, orientation, room reverberation, wall reflections, room structure etc. [...] RIR is analogous to AN IMPULSE RESPONSE in LTI systems. The IR is considered to characterise the behaviour of an LTI system completely.....". Given that you are using FIR filters everywhere in the paper...Shouldn't this part be rephrased / removed completely?

The 350 tap LMS is FIR
RIRs are FIR.
Equations (1) and (2) are FIR

THE NATURAL ROOM FILTER. A room is already a filter. When you pump sound in the room, as soon as you have reflections, you get interference. Interference puts a grid of "nodal points" in the room. There are certain locations in the room that if you put a microphone there, naturally, you will be recording a very low signal. This is the reason behind positioning speakers at certain distances from the walls so that you have the waves projected at the back of the speaker work to your advantage. I am looking at your Fig 9 and I am thinking that I can almost see the grid in the Y direction. The point is: Given a room, a mic and a speaker, the mic might be recording reduced sound levels because of its natural positioning. WITHOUT ANC. The room itself is filtering. Fig 9's caption says "...It can [BE] seen that the good anti-noise source positions..."...perhaps these are not "good", but "natural" (?). Because of the physics...That's where you expect them to be. Not because of the ANC. So, what you are saying is, I am taking into account the natural acoustics of the room to correct for what the signal sounds like at the point of the mic knowing the room it propagates and this helps in ANC (?)...
Fig 11 shows that the ANC is working. This is an ANC, clearly targeting the band $(0 \ldots \approx 1000Hz]$. What is its effect on useful signal that occupies the same region?
SIMULATION METHOD / NOVELTY:. Figure 3 implies that you are using the ray tracing model of sound propagation, which, for a rectangular empty enclosed space is fine. But we are not in the 70s. So, in the presence of this fantastic (and FOSS) piece of software, maybe you could use impulse responses from realistic (and large) auditoriums with specific constraints (?). In i-Simpa you can control everything...And it's free...Which means that you can download it now and start playing.

Hope this helps.


In most implementations the FFT returns the following for the DFT:
$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j 2\pi nk/N}$$
Which would result in the total output being scaled by $N$.
The FFT implementation may or may not be scaled by N so it is always good to check the documentation directly for the fft that is being used. The total power is the sum square of all the FFT bins for a full FFT of length N where the time series data is of length N, once proper scaling is accounted for an assuming no windowing beyond the implied rectangular window is used. (I refer to another post that covers how to properly handle the windowing for power estimation).  So for the result as given above, the total signal power would be (and assuming units of voltage or current in the time domain, with a normalized resistance of 1 then the units of the output would be watts):
$$P = \sum_k\bigg(\frac{X[k]}{N}\bigg)^2$$
Besides reading the documentation, a quick test is to do fft([1, 1, 1, 1]) where fft() refers to whatever fft function is used; if the function needs to be scaled by $N$ before computing the total power, the result for above would be [4,0,0,0]. If the result was [1,0,0,0] then the DFT would have already been scaled by $N$ or another option is the FFT is scaled by $1/\sqrt{N}$ so that the fft and ifft are symmetric operations and the result would be [2,0,0,0].
Tne numpy fft function the OP is using rfft only returns the positive half spectrum, which would be the complex conjugate of the negative half spectrum for real signals and handles odd and even $N$ differently in terms of how many samples are included. For a consistent result, I recommend using fft which will return $N$ samples scaled by $N$, in which case to estimate the power you would first divide the fft results by $N$ and then sum the square of all the bins (consistent with Parseval's Theorem).
Also for Python it is recommended to use the fft in scipy.fftpack and not numpy (which is now in both for reverse compatibility reasons but is expected to eventually be removed from numpy; numpy is to be used for the lower level vector processing and broadcasting operations that scipy depends on but not the higher level numerical processing algorithms that are now all in scipy).
If windowing is applied, then the power computation needs to consider the coherent and non-coherent gain of the window. This is straightforward for a single tone whose power spectral density only occupies one bin; it will be reduced by the coherent gain of the window which is given by the sum of the window coefficients and scaled by $N$ (which is always a loss). However for signals with spectral content that spans multiple bins, we must be careful to consider that the bandwidth of a windowed DFT bin will always span multiple bins, so after the decrease due to coherent gain, there will be an increase due to the double counting as all the bins are summed in the power computation previously given. The net result is a processing gain that is applicable to signals with spectral densities spanning multiple bins and described further in the links below.
Find the Equivalent Noise Bandwidth
Effect of windowing on noise


Moving average filter using DF-I and DF-II

To start, notice that the diagrams are in the $\mathcal{Z}$ domain, multiplying by $z^{-1}$ is the delay of one sample.
DF-I
Writing the equation for $y_n$ of the DF-I
$$y_n = \sum_{k=0}^{M} b_k x_{n-k} - \sum_{l=1}^{N} a_{l} y_{n-l}$$
Comparing with the equation of the moving average
$$y_n = \sum_{k=0}^{M} \frac{1}{M+1} x_{n-k}$$
we see that we could realise the moving average by settings $b_k=1/(M+1)$, and $N$ set to zero.
DF-II
For the DF-II realization you the sequences $a$ and $b$ must have the same size, so what we can do is to set $a_l=0$, and $b_k = 1/(M+1)$.
If you choose to compute a moving sum then $b_k=1$ and you can save the multiplications, the moving average still can be computed by a gain, thus, performing one multiplication per sample instead of $M+1$ multiplications per sample.


You can increase the sampling rate, instead of the signal $x(nT_s) = x(n/F_s)$, you can play it at say $\frac{11}{10}F_s, \frac{12}{10}F_s, \text{or}\ 2F_s$. So you have:
\begin{align}
x\left(\frac{n}{\frac{11}{10}F_s}\right) &=x\left(\frac{10n}{11F_s}\right) = x\left(n\frac{10}{11}T_s\right)\\
x\left(\frac{n}{\frac{12}{10}F_s}\right) &=x\left(\frac{10 n}{12F_s}\right) = x\left(n\frac{10}{12}T_s\right)\\
x\left(\frac{n}{2F_s}\right) &=x\left(\frac{1}{2}\frac{n}{F_s}\right) = x\left(n\frac{1}{2}T_s\right)
\end{align}
In MATLAB for instance, you can try the code below.

clear all

load handel.mat
sound(y, Fs) % original
sound(y, 1.1*Fs) % at 1.1*Fs
sound(y, 1.2*Fs) % at 1.2*Fs
sound(y, 2*Fs) % at twice Fs

Gradually increasing the frequency to hear the changes in speed. You can try even higher, like 4 times, to hear the differences. The handel.mat file should be in there by default.

EDIT:

You $x(t)$ in your program is samples of the signal for a certain time duration, check how many samples you have per second (i.e. the sampling rate or frequency $F_s$), and $T_s = 1/F_s$. And $x(t) = x(nT_s) = x(n/F_s)$ with $n = 0, \ldots, N-1$ where $N$ is the number of samples. Taking a higher $F_s$ corresponds to taking smaller $T_s$.


You're probably better off skipping the $5u(t)$ idea and going straight to a periodic input, probably a sinusoid.  A phase locked loop is for signals that have a phase to lock onto, whereas $u(t)$ doesn't have a phase.  The Stanford reference looks pretty good.  Start by comparing that to your feedback system:  Note that the phase detector is a multiplication, where you appear to have an addition.  Also notice that the VCO isn't actually putting out $\hat \theta$, but $\hat s(t) = \cos\left(\omega_{lo}t + \hat\theta(t)\right)$.

So after the multiplication you end up with the sum of two new sinusoids, one at the sum of the frequencies of $\hat s(t)$ and $s(t)$, the other at their difference.  The loop filter is there to remove the higher frequency one, leaving just the difference, which will be $\sin\left(\theta(t) - \hat\theta(t)\right)$.

Now consider a few different situations.  First, the easy one - the input is running at exactly the same frequency as the local oscillator, so $s(t) = \cos(\omega_{lo}t)$.  Assume there is nothing on the input to the VCO to begin with, and so $\hat\theta(t) = 0$ too.  The output from the phase detector will be only the high frequency term, and that will get filtered out by the loop filter, leaving $e(t) = 0$.  The system is in equilibrium, everyone can go home and have a beer.

But now the frequency of the input signal increases a little.  $\theta(t)$ is now a linearly increasing value (if it were a constant, then the input signal would still have the same frequency, but just be delayed or advanced a bit).  Now the phase detector is also producing a low frequency term, $\sin(\theta(t))$ in addition to the high frequency one.  This is passed through the filter, and since the frequency only increased a little, $\theta(t)$ is small so $\sin(\theta(t)) \approx \theta(t)$.

This means that $e(t)$ is no longer zero, and so the phase of the VCO begins to increase; $\hat s(t)$ runs faster.  If it hasn't caught up to the frequency of $s(t)$, then the phase difference will continue to increase, and since that difference is being integrated by the VCO, the VCO frequency will increase faster.  Eventually the VCO will be running at the same frequency as $s(t)$.  However, note that it will be lagging in phase:  To run at the same frequency, $\hat\theta(t)$ needs to be increasing at the same rate as $\theta(t)$.  For $\hat\theta(t)$ to be increasing linearly, $e(t)$ needs to be a constant, and non-zero, which means $\theta(t) - \hat\theta(t)$ must be non-zero.

If the frequency of the input signal were to decrease, you would get a similar situation, except that now the VCO output would be leading the input signal rather than lagging.

In a real situation, the input signal frequency is usually changing all the time (otherwise you wouldn't need a PLL to track it), and all this is happening constantly.  Hopefully you can see from this description that if the input frequency changes too quickly then the system can be thrown off (remember the small angle approximation I used).

I believe the principles for a discrete-time PLL are similar.  Sampling rates with respect to the input frequency will need to be considered, as usual.

Hope that helps, I'm not sure if you are completely familiar with continuous PLLs and just looking for the differences when you go discrete, or the whole thing is a little fuzzy for you.  Since you've got continuous-time looking variables in your feedback system diagram I've assumed the whole PLL game is of interest to you.


They say that $n(t)$ is white noise with a double-sided power spectral density (PSD) of $N_0/2$, i.e., the PSD is given by
$$S_n(f)=\frac{N_0}{2}\tag{1}$$
The auto-correlation function is the inverse Fourier transform of $(1)$, which is
$$R_n(\tau)=E\big\{n(t)n(t+\tau)\big\}=\mathcal{F}^{-1}\big\{S_n(f)\big\}=\frac{N_0}{2}\delta(\tau)\tag{2}$$


Coordinate system in digital image processing

Different operating systems, standards, and graphical softwares implement distinct (but quite related) conventions for representing, storing, manipulating, or displaying two-dimensional data on computers. The distinction is about the orientation of the axes.
Basically we have 4 main object types :

$f(x,y)$ , x-y continuous, is a function of mathematical relevance, such as an analog image.
$f[n,m]$, n-m integers, is a sequence obtained by sampling $f(x,y)$, or  synth-generated.
$F(i,j)~$, i-j integers, is a matrix to represent the data in its rows-i and columns-j.
$f[i,j]~~$ , i-j integers, is an array to store the data used in a computer program.

The arguments $x,y,n,m,i,j$ can also be replaced with $x_1,x_2,n_1,n_2,i_1,i_2$, the subscripts indicating the dimension order of the variables, such as the first or the second dimension.
A sequence $f[n_1,n_2]$ can represent a 2D digital image, or a mathematical function defined over the coordinate system of axes $n_1,n_2$. Its samples are placed in the correct elements of the associated matrix $f(i,j)$ for the operation to produce the desired effect at the output.
MATLAB fundamentally uses MATRIX based dimension ordering. The first sample A(1,1) is at the top-left corner, first dimenison $i$ is along vertical-down the rows, and the second dimension, $j$, is along the columns (horizontal-right).
This is equivalent to placing a coordinate system of $n_1 ~, ~n_2$ with its origin (0,0) at the element A(1,1); the first axis $n_1$ points along the first dimension $i$ (vertical down); and the second axis $n_2$ points along the second dimension $j$ (horizontal right) of the matrix A(i,j).
This mapping is obtained when you rotate the conventional coordinate system ($n_1$ horizontal-right, and $n_2$ vertical-up, origin at bottom-left) by 90 degrees clockwise with respect to its origin (0,0) and placing the origin at the element A(1,1) of the matrix.
In a parallel fashion, the 2D-DFT (discrete Fourier transform) $F[k_1,k_2]$ of the sequence $f[n_1,n_2]$ is computed by the statement:  F = fft2( f ) which places the first dimension, $k_1$, along the rows $i$, and the second dimension, $k_2$, along the columns $j$ of the output matrix $F(i,j)$ to represent the DFT sequence $F[k_1,k_2]$. Which is also alligned with the first and second dimenions of the input matrix $f(i,j)$ that represents the seqeunce $f[n_1,n_2]$.
Below is an oldskool discussion of a few mapping modes, and functions in MATLAB related with data orientation. As long as mappings are used consistently, they all yield the same results when interpreted correctly.
% SEQUENCES, BMP IMAGES, MATRICES and MATLAB FUNCTIONS :
% ------------------------------------------------------
%
% All data processing and display in Matlab is done via MATRICES A(i,j).
% But, theory of image processing is based on SEQUENCES f[n1,n2], F[k1,k2].
% When processing images, orientation of axes become relevant across
% functions such as CONV2(), FFT2(), IMSHOW(), STEM3(), SURF().
% And a MAPPING convention from f[n1,n2] into A(i,j) should be used.
%
%
% A matrix A(i,j) is indexed by vert rows i, and horz columns j.
% A sequence f[n1,n2] is indexed by horz-right n1, and vert-up n2.
% Which is the most typical, and natural, orientations for n1 and n2.
%
%
% The mapping convention depends on the functions being called :
%
%           1- PROCESS functions: conv2(), fft2(), filter2()
%           2- DISPLAY functions: imshow(), surf(), stem3()
%
%
% We consider following mapping modes between the samples of a sequence
% f[n1,n2] and elements of a matrix A(i,j) as follows:
%
%
% MM-0 : 90d CW ROTATED MAPPING :
% ----------------------------------------------------------------------
%  f[0,0] --> A(1,1), and "n1" grows DOWN from the top-row of A.
%
%    n1 = i-1   , n2 = j-1     ===>          A(i,j) = f[i-1, j-1]
%    i  = n1+1  , j  = n2+1    ===>        f[n1,n2] = A(n1+1,n2+1)
%
%  ---o---1--------2-----> j (n2)
%   1 | f[0,0]  f[1,0]
%   2 | f[1,0]  f[1,1]    Ex mapping of f[n1,n2] into 3x2 matrix A(i,j)
%   3 | f[2,0]  f[2,1]
%     |
%   i v
%  n1
%
%
%
% MM-1 : FLIP-DOWN MAPPING :
% ----------------------------------------------------------------------
% f[0,0] ---> A(1,1) and "n2" grows DOWN from the top-row of A.
%
%    n1 = j-1   , n2 = i-1     ===>          A(i,j) = f[j-1, i-1]
%    i  = n2+1  , j  = n1+1    ===>        f[n1,n2] = A(n2+1,n1+1)
%
%  ---o---1--------2-----3---> j (n1)
%   1 | f[0,0]  f[1,0]  f[2,0]
%   2 | f[0,1]  f[1,1]  f[2,1] Ex mapping of f[n1,n2] into 3x2 matrix A(i,j)
%     |
%   i v
%  n2
%
%
% MM-2 : BMP MAPPING :
% ----------------------------------------------------------------------
% f[0,0] is stored in A(N2,1) and "n2" grows UP from the bottom-row of A
%
%    n1 = j-1   , n2 = N2-i     ===>          A(i,j) = f[j-1, N2-i]
%    i  = N2-n2 , j  = n1+1     ===>        f[n1,n2] = A(N2-n2,n1+1)
%
%  n2 ^
%   2 | A(1,1)  A(1,2)         An example mapping into 3x2 matrix A
%   1 | A(2,1)  A(2,2)
%   0 | A(3,1)  A(3,2)
%  ---o---0--------1-----> n1 (j)
%
%
%
% Axis orientations of functions CONV(),FILTER(),FFT2():
% ----------------------------------------------------------------------
% They assume 90d ROTATED / or FLIP DOWN index mapping by default.
%
%      A(1,1)---------> n2,k2,j (horizontal - column variable)
%           |
%           |  A(i,j)
%           |
%           v n1,k1,i (vertical - row variable)
%
%
% NOTE: Strictly speaking, MATLAB does not care which mapping was used to
% generate the matrix A(i,j). Rather it treats "i" as the first dimension
% and "j" as the second dimension, and outputs acordingly.
%
%
% Axis orientation of DISPLAY functions STEM3(),  SURF()            :
% -------------------------------------------------------------------
% It produces the plots according to NATURAL X-Y orientation in which
% n1 point horizontal right, and n2 point verticcal up and (0,0) being
% at the bottom. But this requires that the MATRIX A(i,j) was filled in
% according to FLIP-DOWN mapping mode.
%
%                          A(N2,N1)
%                          f22
%                       f12   f21
%     A(N2,1) (n2-i) f02  f11  f20 (n1-j) A(1,N1)
%                      f01  f10
%            n2-axis      f00       n1-axis
%                        A(1,1)
%
%
% Image display function IMSOW(A(i,j)) assumes a BMP mapping:
% -------------------------------- --------------------------
% IMSHOW displays the matrix contents in its row-column order.
% Therefore if you want to display a sequence f[n1,n2], or F[k1,k2] using
% the IMSHOW() function, then in order to have the correct orientation
% according to natural x-y coordinates, you should use BMP mapping mode
% on the matrix which is to be displayed.
%

SUMMARY
If IMSHOW() iwill be used to display images, or their FFT results, then correct alignment with $n_1,n_2$ and $k_1,k_2$ requires BMP based mapping to be used to fill in the associated matrix. If STEM3() or SURF() will be used to get a 3D display of the sequences, or FFT results, then flipp-down mapping produces correct orientation. Other functions CONV2(), FFT2() etc., work equally well with either Rotated or Flip-Down mapping modes.


Designing digital low pass filter with low pass-band group delay

Minimum phase filters will not give you a near constant group delay. You can design a non-linear phase FIR filter with a linear desired passband phase with a specified group delay that is smaller than the group delay of the corresponding linear phase filter. If you use a least-squares criterion, this is equivalent to solving a system of linear equations.
As an example, take an FIR lowpass with $61$ taps. The group delay of a linear-phase filter would be $30$ samples. Now we specify a group delay of $20$ samples instead:

w = pi * [linspace(0,.2,100), linspace(.25,1,375)];   % frequency vector
D = [exp(-1i*w(1:100)*20), zeros(1,375)];   % desired complex frequency response
W = ones(1,475);   % unity weighting
h = lslevin(61,w,D,W);   % least squares design

The function lslevin.m can be found here.
The result looks like this:

Note that the passband group delay is of course not exactly constant, but the desired group delay of $20$ samples is approximated quite accurately over a large part of the passband.
Trying to achieve a low group delay in the passband comes with a penalty in the magnitude response. The figure below shows the magnitude and group delay responses of $4$ length $61$ FIR lowpass filters with different desired group delays in the passband. One is a linear phase filter with a constant group delay of $30$ samples, and the other $3$ filters are non-linear phase low delay filters with specified passband group delays of $25$, $20$, and $15$ samples respectively. Clearly, the magnitude response deteriorates with decreasing desired passband group delay.


It's pretty simple if you know a few basic things:

a single pole on the real axis corresponds to an exponential sequence. The elements of that sequence either have the same sign, or their signs alternate, depending on the angle of the pole ($0$ or $\pi$). This should allow you to pair figures $(c)$ and $(e)$ with figures $(1)$ and $(6)$.
a complex conjugate pole pair corresponds to a sinusoid with exponentially decaying or increasing amplitude, depending on whether the poles are inside or outside the unit circle. This should help you pair figure $(f)$ with the corresponding sequence.
if there are only zeros in the PZ-diagram (and poles at the origin or at infinity), then the corresponding sequence has finite length. Note that since all shown sequences are causal, there must be the same number of poles at the origin as there are zeros. Those poles are not shown in the figures, which could be considered bad practice.
linear phase sequences have certain symmetry properties in the time domain as well as in their PZ-diagrams.
the formula for the geometric series helps with figure $(5)$.

Also take a look at this related question and its answers.

  a5,b2,c6,d4,e1,f3


I would search for "segmentation benchmark dataset" or so. One of the first results seems to be useful for your task:

https://www.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/segbench/BSDS300/html/dataset/images.html


Is it possible to use k-means algorithm with just one vector

Certainly possible. You can think of it in the way you described but I would say it just becomes an adaptive thresholding assuming a bimodal distribution: We admit that there are two main modes in the data and we try to identify those by grouping the values around the means of the clusters. These groups form clusters. The means are positioned such that each element contributes to the average of the cluster it belongs and thus creating a mode nearby.
It is almost sure that the minimum and maximum will go to separate clusters.


Model Validation After Estimation for System Identification Task (Assistance with MATLAB Code)

I agree with AlexTP's idea of simply plotting the instantaneous error between the model's output and the system's output over time. The exact implementation is a bit tricky though, because you are not only estimating the impulse response (using LMS) but also doing an "inverse" filter operation to estimate the unknown inputs. One way is to use a standard machine learning approach outlined in this answer.

Split your known input data $\mathbf d$ into two halves - training and testing (or validation$^*$).
Using the training half, run LMS and estimate the channel coefficients $\hat{\mathbf{h}}$.
Now use your test set and apply them to the noisy system to obtain outputs $x_n = \mathbf h^T \mathbf d + e_n$ and also apply them to the filter learned in step 2 to get $\hat x_n = \mathbf{\hat h}^T\mathbf d$. Compute the error$^\dagger$ between $\hat x_n$ and $x_n$ and make sure it's acceptable/small.$^\&$
If you are happy with the performance in step 3, your filter $\hat{\mathbf h}$ is ready for real world where you use the inverse of the filter $\hat{\mathbf h}$ to estimate the unknown inputs.$^{\ddagger}$

$^*$ Since you aren't doing model selection, you don't need a validation set, or in other words, your test and validation set are the same. If you do care about model selection (eg. choosing the best filter order) you can split your data into 3 pieces - training, testing and validation.

$^\dagger$ I don't think there's a single best way of computing this error. You could simply plot $\hat x_n - x_n$ vs $n$ and make sure it stays close to zero. Or perhaps you could aggregate them as a mean absolute error or mean squared error by summing over all $n$.

$^\&$ Alternatively, you can use the inverse filter of $\hat{\mathbf h}$ and estimate the inputs and compare them with the known inputs. Again, there's no single best way of showing error and would depend on your application.

$^{\ddagger}$ Since this phase has unknown inputs you can't really compute a performance metric for your model. Or as this answer puts it "you can only speculate about the quality of your model output using the results of your validation phase".


Deriving Frequency Response for 2-pole Zero-Delay Feedback State Variable Filter

I do not see any reason why this could not be derived algebraically. But this would be a good application for Mason's gain formula; the references below may be helpful for that approach:

Rick Lyon's Tutorial

Mason's Gain Formula - Wikipedia

Youtube tutorial

Tutorial's Point tutorial


So, a sound signal $s$, in the DSP sense, is but a sequence of real-valued digital numbers:

$$ s = (1, 232, -12, 121, 0, 131, \dots) $$

Integrating that is nothing but adding these up, so the integrated signal $r$ would have the formula:

$$ r[n] = \sum\limits_{i=0}^n s[n] $$

or, writing that recursively, is the signal that it was at the last sampling time, plus the next input sample:

$$ r[n] = r[n-1] + s[n]\quad \text{ for }r>0\\
r[0] = 0$$

Now, that is what DSP people call a recursive system, or in this specific case, a single-pole IIR filter (IIR: infinite impulse response, because you send one pulse in, you see it in the output forever). It's transfer function $H(z)$ is

$$H(z) = \frac{1}{1+a\cdot z^{-1}} $$

Your $a=1$, here, because the last value of $r$ doesn't get multiplied with a dampening factor. You will find a lot of literature out there explaining the properties of IIR filters, but for all practical things, you could build a low-pass filter out of your integrator, if you use a factor $a<1$.


How do I attenuate a tidal signal in depth time series data

The answer is "No, you cannot attenuate the impact of tide on pressure derived depth data you have already captured without extraneous information".

The reason for this is that the hydrostatic pressure derived depth measurements from the tracked fish ($depth_{fish}(n)$), is only a component of the captured signal.

The illustrative example here is a "freakishly intelligent fish" which manages to stay at one point (with respect to the bottom of the sea) and adjusts its buoyancy with respect to pressure. When the tidal potential piles up the water in a water column above the fish, it ascends. When the water recedes, the fish descends. The pressure derived depth of this fish is constant, when clearly, it has been going up and down.

The hydrostatic pressure signal on the fish body is $P_{recorded}(n) = P_{tidal}(n) \cdot P_{fishDepth}(n)$ where $n$ is some variable that indicates time.

You are FFTing the $P_{recorded}$, aiming at suppressing the $P_{tidal}$, but, multiplication in the time domain is convolution in the frequency domain. Therefore, just by having $P_{recorded}$, you have both $P_{tidal}, P_{fishDepth}$ intertwined. That's not entirely bad news but it certainly is not the best news ever.

You need to somehow fix one to recover the other and to an extent, whether with freakish frightening accuracy as above or simply casual preference, fish do regulate their buoyancy, but perhaps faster than the speed by which the tidal potential moves over the globe. Therefore, your task might become as easy as applying a low pass filter to attenuate the tidal component which is more likely to exist at the low frequencies. The speed at which the fish are moving versus tidal variation is probably your best "weapon" here. For a bit more background information, you might want to see this link which is not exactly what is asked here but is relevant. (Note: It depends very much on the scale of observation too. Fish might have specific life patterns that are modulated by the time of day. So again, their signal might "hide" in the tide because it just so happens to be driven by the celestial objects that also drive the tide.)

Now, the reason I am asking you about your accuracy specifications is because you might be able to construct some kind of model by which you could isolate the movement of the fish from tide and/or depth.

The simplest way to think about this is to say, "My sensor says 23 meters. Where is the fish? It's at XYZ location. When is the fish? It's at some hh:mm:ss time. Right, what is the depth (from the datum) at the where? It's 50 meters. What is the depth at the where and the when? (Due to tide), it's 52 meters.", Ah! So, the tide has added two meters, therefore the fish is now at 21 meters from the bottom.

It might be possible to work out a probabilistic model too. That answers the question of what is the probability of a depth value, given the prior formed by the sum of bathymetry and tidal data and the data obtained by the sensor.

In any case, you either need extraneous information (bathymetry,tide) or to make an assumption (tidal signal much slower, therefore at very low frequencies, apply a low pass, done).

Finally, please note that bathymetry does play a significant role. For a trivial example, imagine a basin that fills up with water (and fish) during the tide and then preserves water (and fish) even after the tide has receded. For those fish, the readings will always be (for instance) between 4-5 meters when around them, the depth could well be zero. So, your $P_{tidal}$ is actually a sum of bathymetry and tide.

I suspect that this is probably "bad news" (?) but I hope it helps anyway.


That paper is very poorly written. It's no wonder you have a hard time understanding it; I certainly did!

What I think is happening is that several of the DCT coefficients are being thrown out (windowed). So if the window is 100%, then all the DCT coefficients are kept; if the window is 30% then 70% of the coefficients are thrown out (zeroed).

See Figure 11 of that paper:

The left hand signal are the DCT coefficients. The right-hand graph is the applied window. The two are multiplied together.


What is meant by Two pass FIR Filter?? (A basic question)

You can't in general categorize a two-pass filter in any of those listed filter types. "Two-pass" means that the filter processes the data in two passes; the output of the first pass is used as the input of the second pass.

The two-pass filter's impulse response is the convolution of the impulse responses of the filters in the two passes. This comes from the algebraic properties of convolution:

$$\text{output}= f*(f*\text{input}) = (f*f)*\text{input},$$

where $*$ denotes convolution and $f$ is the impulse response of a single pass. The two passes have an identical impulse response to match your description. Because convolution is equivalent to multiplication in the frequency domain, the two-pass filter's frequency response will be the square of the frequency response of a single pass.

With infinite impulse response (IIR) filters, also a configuration where the second filter processes the data backwards can be useful, and the two filters may be connected in parallel so that the outputs from the two passes are summed to form the composite filter output. Such configurations are not useful with FIR filters.


What information does fourier transform carry?

There are a variety of Fourier Transforms (and Series) such as Continuous Time FT, Discrete Time FT, Discrete FT all of which are generally attributed to the fundamental assertion made by J.B.Fourier about at the beginning of 19th century, which claims (without proof) that "if a continuous-time signal (function) x(t) is periodic with T, then it is possible to represent that signal x(t) as an infinite sum of harmonically related trigonometric functions (sines and cosines) as $ x(t)= \sum{ a_k \sin { 2\pi kt \over T} + b_k \cos { 2\pi kt \over T}}$  in which the weights $a_k$ and $b_k$ (the Fourier Coefficients) represent the amount of that particular harmonic in the signal x(t) being analysed"

In fact the above argumentation is strictly for Continuous-Time Fourier Series. But the core idea is generalised into Fourier-Transforms of aperiodic and periodic (with the help of impulse $\delta (t)$ functions) signals. There are conditions on which signals can have such a representation.

In essence, computing a Fourier Transform means finding those coefficients $a_k$ and $b_k$ for which the method is suggested by the analysis equation of the Fourier transform, while the equality in the first paragraph is noted as the synthesis equation.

Eventhough it is quite solid to understand the meaning of those sinusoids inside a periodic signal, when it comes to non-periodic signals, for which we use the Fourier Transforms, the exact meaning of what a single sine wave represents inside such a signal is a little vague and instead we emphasize the transient character of that signal under concern and the necessity of existance of a continuum of infinetely many sine waves.


Because of the symmetry of the $\operatorname{sinc}(t)$ we have
$$R_{hh}(\tau)=\int_{-\infty}^{\infty}\operatorname{sinc}(t) \ \operatorname{sinc}(t-\tau) \ \mathrm{d}t=\int_{-\infty}^{\infty}\operatorname{sinc}(t) \ \operatorname{sinc}(\tau-t) \ \mathrm{d}t = \operatorname{sinc}(\tau) \ast \operatorname{sinc}(\tau) $$

Analyzing in the frequency domain we have
$$\mathcal{F}(R_{hh}(\tau)) = \mathcal{F}(\operatorname{sinc}(\tau))^2=\operatorname{rect}(f)^2=\operatorname{rect}(f)$$

$$\Rightarrow R_{hh}(\tau)=\mathcal{F}^{-1}(\mathcal{F}(R_{hh}(\tau))) = \mathcal{F}^{-1}(\operatorname{rect}(f))=\operatorname{sinc}(\tau)$$

where $\operatorname{rect}(f)$ denotes the rectangular function in the frequency domain


Why can't realisable IIR filters have linear phase?

The frequency response of a real-valued discrete-time system with linear phase has the form

$$H(e^{j\omega})=A(\omega)e^{-j\omega\tau},\qquad\omega\in [-\pi,\pi]\tag{1}$$

where $A(\omega)$ is either a real-valued even function or a purely imaginary odd function, and $\tau$ is some real-valued parameter (the delay). If $A(\omega)$ is purely imaginary, then there is a phase jump of $\pi$ at $\omega=0$, but otherwise the phase is linear. This latter type of phase is usually referred to as "generalized linear phase", and it is the type of phase exhibited by filters with odd symmetry in the time domain.

The impulse response of the corresponding system is given by

$$h[n]=(a\star d)[n]\tag{2}$$

where $a[n]$ is the inverse discrete-time Fourier transform (IDTFT) of $A(\omega)$, $d[n]$ is the IDTFT of $e^{-j\omega\tau}$, and $\star$ denotes convolution.

The sequence $a[n]$ is either symmetrical or anti-symmetrical due to the properties of $A(\omega)$ (real or purely imaginary), and hence non-causal. The sequence $d[n]$ is given by

$$d[n]=\frac{\sin[\pi(n-\tau)]}{\pi(n-\tau)}\tag{3}$$

which simplifies to $\delta[n-\tau]$ if $\tau$ is an integer.

Consequently, if $a[n]$ has infinite length (as is the case for IIR filters), $h[n]$ cannot be causal and have linear phase at the same time. Only if $a[n]$ has finite length (thus FIR) can it be made causal by shifting to the right. That's why FIR filters can have linear phase whereas causal IIR filters cannot.


Image processing for extremely small feature detection

I'm going to go ahead and assume currently using Matlab implies you're on a laptop or PC-alike device with a modern x86 CPU; optimal/efficient/power-saving DSP depends on the hardware you're on, but if we have the features of a fully fledged PC platform, it becomes pretty easy :)

I think this is the classical case where a linear search doesn't hurt.

Assume your pictures come in a in a row-major format, i.e. each pixel lies in memory next to its left and right neighbours, and the top and bottom neighbours are one row width of pixels away.

Since memory bandwidth is the limiting factor for simple operations on data, it's always important to minimize the times you fetch new memory – things that are in the CPU registers can directly be used for calculations, whereas things that still are in RAM easily take a couple hundred multiplications worth of CPU time to just get loaded, which is why there is an intermediate RAM cache in your CPU, which your processor can use efficiently – unless you jump wildly through memory as you process, because then it newer knows what it should keep in cache and what to throw out to make space for new stuff.

So, row-major greyscale.

You go to the first row in your picture, and search for the brightest pixel; you keep a variable storing both the value (necessary for searching) and the index of that pixel (column). At the same time, you keep a sum of pixels. At the end of the row, you store that sum and the index of the brightest pixel.

Move on to the next row, and so on.

At some point, you might want to terminate; either you've run out of rows, or your algorithm has produced enough high row sums, and isn't producing them anymore, so you've probably found your circular spot.

With the brightest pixel column indices and knowing in which rows that happened, you have your region of interest in linear time, without thrashing the CPU cache.

If this isn't fast enough for you, the nice thing about the operations you've been doing (sum, comparison with value) is that modern x86 CPUs have instructions that let you these operations with up to 16 operands at once (depending on the data type you're using underneath), so this might be a good time to look into that. I'd recommend not writing SIMD intrinsics from hand for a start, but looking into VOLK (Vector Optimized Library of Kernels) which should have the routines you need, at least for floating point, already.

Yet another thing: this can be highly scaled on multiple CPUs. If you've got four CPU cores, the first could do the upper fourth of the image, the second the next fourth and so on. Same goes for doing this with OpenCL on your graphics card – but do remember that the whole data marshalling in and out of graphic cards and GPU drivers might eat its performance children rather easily.

So, with your Region Of Interest detected, you can do multiple things. I'd just go ahead and copy the smaller subimage of which I'm sure contains the whole spot into a new, compact image, and probably store it for later verification (at least as long I'm still developing); then I'd go and let one of the classical circle-finding approaches loose on the image: Typically, the circular Hough transform will do just fine. Notice that you could have done that directly with the original image, but you'd have lost a lot of time spent on transforming things that definitely don't contain circles. The OpenCV has circular Hough built-in.

A bit of development advice: don't prematurely optimize! Get your raw ROI detection working with whatever programming language you prefer (since we're talking about most efficient implementations here, this is pretty likely to be C/C++, but needn't be! A functional programming language might make it a lot easier for the compiler and runtime to parallelize things, though it will be harder (if possible at all) to later on sneak SIMD intrinsics into the code.). You'll probably notice that pretty much everything beats the speed of the Matlab interpreter by orders of magnitude, whereas the natively implemented low-level algorithms within Matlab are highly optimized (you pretty much will have a hard time beating the speed of Matlab's FFT implementation – which is FFTw, by the way, or it's matrix/vector operations, which are BLAS/LAPACK, and so on). But practical algorithms not only do elementary operations but also have program flow defined in matlab M-code, and that, performance-wise and from the point of elegance, is stuck somewhere in the early 80's, at times. So you'd be surprised how much you often can win purely by practically writing the same program in Python, using numpy, scipy, pyOpenCV, without touching low-level coding at all.


The first definition works for deterministic as well as for random signals. For random signals we define the autocorrelation by

$$R_x(\tau)=E\{x^*(t)x(t+\tau)\}\tag{1}$$

where $E\{\cdot\}$ is the expectation operator. For deterministic power signals (i.e., signals with finite but non-zero power and infinite energy) we define

$$R_x(\tau)=\lim_{T\to\infty}\frac{1}{2T}\int^{T}_{-T}x^*(t)x(t+\tau)dt\tag{2}$$

From $(2)$ it is clear that for deterministic power signals both definitions of average power are equivalent, because for random signals as well as for deterministic power signals we have

$$S_x(f)=\mathcal{F}\{R_x(\tau)\}\tag{3}$$

Also take a look at this related answer.

You are right about the power of white noise. In continuous time, white noise has infinite power. Only discrete-time white noise has finite power. In sum, teachers and professors are only humans.


How to recover a signal from its complex envelope

According to the usual definition, if you have a complex low pass signal (complex enevelope) $u(t)$ and a reference frequency $f_c$, then the corresponding pass band signal is

$$u_p(t)=\text{Re}\{u(t)e^{j2\pi f_ct}\}\tag{1}$$

So if $u(t)$ is given by

$$u(t)=A(t)e^{j\phi(t)}\tag{2}$$

the real-valued pass band signal would be

$$u_p(t)=\text{Re}\{A(t)e^{j\phi(t)}e^{j2\pi f_ct}\}=A(t)\cos(2\pi f_ct+\phi(t))\tag{3}$$


From what I can glance :

The chattering seems to be high-frequency compared to your signal of interest, it should not be hard to filter this chattering noise. You simply need to identify the frequency band of your signal and the frequency band of this noise.

Could your perform an FFT to analyze the frequencies of your noise? Then design a filter that will not attenuate your signal of interest while attenuating the noise to an acceptable level.

Is it real-time processing or post-processing?


Same z transformed function, but different answers of inverse z transform?

Here's a much simpler method. We have

$$E (z) := \frac{1}{z+4} = z^{-1} \left(\frac{1}{1 - (-4) \, z^{-1}}\right)$$

which we assume to be the $\mathcal Z$ transform of a causal discrete-time signal. Consulting a table of $\mathcal Z$ transforms, we conclude that

$$\mathcal Z^{-1} \left\{ \frac{1}{1 - (-4) \, z^{-1}} \right\} = (-4)^n \, u (n)$$

where $u$ is the discrete-time step function. Hence, the inverse $\mathcal Z$ transform of $E$ is simply the discrete-time signal above after a unit delay, i.e.,

$$e (n) = (-4)^{n-1} \, u (n-1)$$

which is zero at $n=0$. Hence, your 1st candidate answer is (almost) correct. Your 2nd candidate answer is only correct if $n \geq 1$. When working with functions, do not forget what their domains are.


For a system to be causal, number of finite zeros <= number of finite poles. Why?

If the number of finite zeros is not greater than the number of finite poles then the transfer function is proper, i.e., the degree of the numerator polynomial is not greater than the degree of the denominator polynomial.

If the degree of the numerator polynomial were greater than the degree of the denominator polynomial, we would get at least one pole at infinity. This violates the requirement that a causal and stable system must have all its poles in the left half-plane (for continuous-time systems), or inside the unit circle (for discrete-time systems). Consequently, a causal and stable system must have a proper transfer function. This is true in continuous time as well as in discrete time.

These two answers are related: answer 1, answer 2


Why do we need >= 95% throughput for LTE Sensitivity measurement according to 3GPP

It's a definition. There's probably no specific reason
other than "the committee agreed on it".
Point is: this is for comparability of other tests. So, the actual number 95% is less important than that all use the same.