"Does audio signal [processing] specifically have any different filtering method?"

well, there are IIR filters that can be implemented in various forms (like Direct I, Direct II, the transposed Direct I or II, Lattice, "State-Variable", Gold-Rader...), there are FIR filters and FIR filters can be implemented the various ways (transversal or the simple summation, Lattice, and "Fast convolution" using the FFT).

any of these filtering methods used on other signal classes can also be used for audio.


Note that the condition on the relation between the bandwidth of $m(t)$ and the modulation frequency ensures that $h(t)$ is a band pass signal. It can be written as

$$h(t)=\frac12 m(t)e^{j2\pi \nu_c t}+\frac12 m(t)e^{-j2\pi \nu_c t}\tag{1}$$

From our assumption $\nu_m<\nu_c$ it follows that the two expressions on the right-hand side of $(1)$ do not overlap in the frequency domain. Consequently, the analytic signal is simply given by the first term on the right-hand side of $(1)$ (scaled by a factor of $2$), which is the part of the spectrum of $h(t)$ occurring at positive frequencies:

$$h_a(t)=m(t)e^{j2\pi\nu_c t}\tag{2}$$

The Hilbert transform of $h(t)$ is given by the imaginary part of $(2)$:

$$\hat{h}(t)=\text{Im}\{h_a(t)\}=m(t)\sin(2\pi\nu_c t)\tag{3}$$

EDIT: The equation for $h_a(t)$ given in your question is also correct. If you want to continue that way, it's important to see that due to $\nu_c>\nu_m$, $M(\nu+\nu_c)$ equals zero for $\nu>0$. So the integral can be written as

$$h_a(t)=\int_0^{\infty}M(\nu-\nu_c)e^{j2\pi\nu t}d\nu\tag{4}$$

Likewise, $M(\nu-\nu_c)$ is zero for $\nu<0$, so we can extend the lower integration limit in $(4)$ to $-\infty$:

$$h_a(t)=\int_{-\infty}^{\infty}M(\nu-\nu_c)e^{j2\pi\nu t}d\nu=\mathcal{F}^{-1}\{M(\nu-\nu_c)\}=m(t)e^{j2\pi\nu_c t}\tag{5}$$

which is of course identical to $(2)$.


In image denoising far more important then the noise distribution is the noise spatial correlation properties and the prior about the image.

Let's try building some cases and dealing with them.
The model is:

$ y = x + n $

Where $ x $ is the clan image, $ n $ is the Poisson Noise (With mean $ \lambda $) and $ y $ is the noisy image.

Noise Is Poisson Distributed and Spatially White

In this case each pixel $ {y}_{i} $ is distributed with Poisson with mean $ {x}_{i} + \lambda $.

Deriving the the estimator for $ {x}_{i} $ will yield it is $ {y}_{i} $. So no gain here what so ever.

Noise Is Poisson Distributed and Spatially White with Image Prior

If we assume some kind of a prior on the image (Piece Wise Smooth) we can do something.
The simplest idea is use the mean of a neighborhood we believe should have the same value (Think of capturing a wall with the same color all over it).

A very popular and effective way to enforce Piece Wise model is using the Total Variation prior on the Gradients of the image.
Then the MAP becomes something to work with (See Noise, Image Reconstruction with Noise (EE367/CS448I: Computational Imaging and Display, Class 10, Gordon Wetzstein, Stanford University)).

Noise Is Poisson Distributed and Spatially Correlated

The usual technique to deal with spatially correlated noise is working in a Multi Scale manner.

Noise Is Poisson Distributed and Spatially Correlated with Image Prior

Now combine all methods and voila, you have a real world denoiser.


Computer Vision/Image Processing Riddles

Nice idea :)

A few years ago, the Instagram team proposed the Unshredder challenge. The chalélenge is over now, but the game is still interesting and could be extended to more difficult puzzles.


How to take inverse fft of windowed and callibrated fft data?

When using quantized data formats, applying a window function is a lossy process, and arbitrary original data cannot be recovered, either before or after an FFT/IFFT process.

You could try multiplying by the reciprocal of the window (being careful to avoid divide-by-zero), and get a noisy reconstruction (except perhaps for an end point or two).  And noisier at the ends, where the window tapers to near zero, thus potentially adding a huge amount of quantization noise.

Or you could overlap Hamming/Von_Hann windows by 50%, where the sum of the sequence of overlapped windows has a gain of 1.0 (except at the very ends), and thus is not lossy (or as lossy, but still with some amount of tiny numerical or rounding noise).


There are some important clarifications out of this question that are good to point out:

The frequency resolution of a block of data is 1/T where T is the length of the data in time (in seconds). Since the sampling rate and block length are related as

$T = N/f_s$

Where $f_s$ is the sampling rate, it follows that the frequency resolution for a block of N samples will be $1/T = f_s/N$.

Zero padding does not change the frequency resolution, it only interpolates more samples of the Discrete Time Fourier Transform. See these posts:

What happens when N increases in N-point DFT

What proportion of a padded FFT should be actual values

When you do fft(signal,N) in Matlab, N must be larger or equal to the signal length and when larger it will simply append zeros, so identical to zero padding.

So to achieve a 10 Hz resolution bandwidth, you either need to reduce the number of samples, as in fft(sig(1:m)), or reduce the sampling rate, but regardless of sample rate, the length of the block in your fft must be 1/10 = 0.1 seconds long.

Note that windowing if performed will expand the frequency resolution; the finest resolution is achieved with no windowing and is $f_s/N$

Refer to this paper by fred harris which details the resolution bandwidth for various windows: fred harris on the use of windowing


To quote my dissertation:

  A collection of transforms are given the name constant Q and are
  similar to the Fourier transform.

  Computation of the discrete Fourier transform can be very efficient
  when employing the use of the fast Fourier transform. However we
  notice that energy of a signal is divided into uniformly sized
  frequency buckets across the spectrum. While in many cases this is
  useful, we notice situations where this uniform distribution is
  sub-optimal. An important example of such a case is observed with the
  analysis of musical frequencies. In Western music, the frequencies
  that make up the musical scales are geometrically spaced. We therefore
  see that the map between frequency bins of the discrete Fourier
  transform and the frequencies of musical scales is insufficient in the
  sense that the bins match poorly. The constant Q transform addresses
  this issue.

  The aim of the constant Q is to produce a set of logarithmically
  spaced frequency bins in which the width of the frequency bin is a
  product of the previous. As a result we may produce an identical
  number of bins per musical note across the audible spectrum, thus
  maintaining a constant level of accuracy for each musical note. The
  frequency bins become wider towards the higher frequencies and
  narrower towards the lower frequencies. This spread in the accuracy of
  frequency detection closely imitates the manner in which the
  human-auditory system responds to frequencies.

  Additionally, the close matching of notes in western scales renders
  the constant-Q particularly useful in note detection; identifying a
  musical note value rather than an explicit frequency value.
  Furthermore the constant Q simplifies the process of timbre analysis.
  The frequencies of a musical note played by an instrument are often
  comprised of harmonically related partials. The timbre of the
  instrument can be characterised by the ratios of the harmonics. With
  the constant Q transform, the harmonics are equally spaced across the
  bins regardless of the fundamental frequency. This greatly simplifies
  the process of identifying an instrument playing a note anywhere in
  the scale simply by shifting the characterisation across the bins. A
  potential downside to using the constant Q transform is that it
  demands more computation than the Fourier transform.

  An efficient algorithm for transforming a discrete Fourier transform
  (which may be computed with the FFT) into a Constant Q transform is
  detailed in Brown and Puckette (1992).


What is the zero in this transfer function?

$$2s+1=2\left(s+\frac12\right)$$
That's all I can say.


Nyquist Theorem - Why unique frequencies upto Fs/2 and not Fs? f+Fs is start of Aliasing

Since you stated you are still confused, the following graphics may help further understand the other two good answers:

First be sure that you understand what the sampling process (which is multiplying an analog signal by a stream of impulses in the time domain) looks like in both time and frequency; a stream of impulses in the time domain is a stream of impulses in the frequency domain. Multiplication in the time domain is convolution in the frequency domain, so if we sample by multiplying by an impulse every T seconds, we convolve our waveform with a frequency impulse spaced every 1/T Hz:

An example of this is sampling a 3 Hz cosine (where a 3 Hz cosine is represented in the frequency domain by an impulse at + and - 3 Hz) with a sampling clock of 20 Hz as shown in the graphic below. The top spectrum represents the input (analog) spectrum of the cosine wave. The middle spectrum is the spectrum of our sampling process (an impulse every 20 Hz including DC), and the bottom spectrum is the digital spectrum resulting from convolving the previous two.

Note that the analog spectrum was replicated in this process; so everything from $-F_s/2$ to $+F_s/2$, what we call the "first Nyquist zone", is duplicated centered on every multiple of $F_s$. (Where $F_s$ is the sampling rate). Thus the unique spectrum is completely defined by what is between $-F_s/2$ to $+F_s/2$, but we could easily show instead the spectrum for 0 to $F_s$, so it is somewhat arbitrary what slice of the spectrum we choose.

Note for real signals, the negative frequency components are the complex conjugate of the positive frequency components, so in that case not unique: they are completely defined by the positive frequency components so in that case that is all we need to present to provide all relevant information (the spectrum from DC to $F_s/2$).

That is not the case for a complex signal, as explained by hotpaw2 and shown in the graphic below.


First of all, a link to an example of the Goertzel algorithm.

Second, the goertzel algorithm is so fast that you don't need to get into multithreading.

Just create a class that does the goertzel algorithm for a frequency and sampling rate that you define when you create an instance.
Make a list or array of instances of that class, then for every sample, run through that array and have it calculate the goertzel for that sample.

In your case, you just need two instances so you don't really need an array of instances.
Fast and simple.
I've done this very thing with several hundred instances on a fairly old PC, and updating the GUI causes more of a load than the processing.


Computational Complexity of Polyphase Resampling

Rational Resampling

10 kHz -> 300 Hz is a rational resampling with relatively benign factors:

$$ \frac{300\,\text{Hz}}{10\,\text{kHz}}=\frac{3\cdot10^2}{10^4}=\frac3{100}\text,$$

meaning that you'd interpolate by 3 and decimate by 100. That's exactly what you're achieving using resample_poly:

Polyphase Decimators

Remember that for decimation by 100 without aliasing, you need to filter your signal to $\frac1{100}$ of it's original bandwidth, which might take quite a lot of filter taps. Sadly, you throw away 99 of 100 samples that you filter in the following decimation step.

Using polyphase techniques, you can cleverly rearrange the $\frac1{100}$-band filter taps so that you don't even do these 99 calculations – and reduce the effort of your decimation filter by the decimation factor, here, 100.

This basically happens by taking your long filter and "deinterleaving" it to 100 disjunct sets, forming 100 polyphase components of your original filter. If you put in the input samples in all of these partial filters, and sum up the output samples of these component filters correctly, you get exactly the same filter. It's really just a way of "arranging" your filter differently. However, if you know which output samples you're going to throw away, you can omit pushing in each input sample to 99 of 100 polyphase filter components – so that you instead "round-robin" the input samples to each of the polyphase filter components. That saves 99% of calculations!

Polyphase resamplers

But what about the upsampling by 3?
Well, another beautiful aspect of polyphase rational resamplers is that you can show that the order of "interpolate by 3 first, eliminate images by $\frac13$-band filter, then apply $\frac1{100}$-band filter to avoid aliasing, then throw away 99 of 100 samples" can be reordered to decimate-filter-interpolate. That's awesome because:

the long filter runs at the lowest rate in the system
you can omit the shorter filter, because it's far more relaxed than the longer one, anyway.

Why is that faster than the $\mathcal O(N+Nlog(N))$ of an FFT filter?

Aside from the aforementioned dramatic reduction in number of necessary operations, remember that even for the shorter polyphase component filters, you can apply FFT filters, if they are long enough so that this pays; in my experience, for rates as benign as yours, that's usually not the case. Compare a post by Tom Rondeau about comparison of FFT FIR filters to polyphase FIRs.

Generally, the beauty of FFT filters lies in the fact that you can reduce the number of operations you need to do to convolve a piece of signal with a filter tap vector. On the other hand, you buy that algorithmic elegance with totally non-linear memory accesses when computing that FFT. Now, if your filter is long, so must be your FFT, and your CPU can't keep all of both the signal vector and the tap vector in localized caches. Rule of thumb says that a single fetch from RAM is worth 200 to 500 Floating point operations[citation needed] on modern x86 computers.

However, even if the amount of operations you need to handle is quite a lot larger in the non-FFT FIR filter, your memory access is (if your FIR implementation isn't totally braindead) inherently linear – and your CPU and its memory controller are really good at fetching linear memory, and even pre-fetching it. So, that works out in your favor quite a long time, until the complexity advantage of the FFT filter really beats the crap out of the naive FIR.

Another aspect of that linearity is easier vectorization using SIMD instructions – if you want to feel disappointed by scipy, install GNU Radio (easy if you're on current Linux distros or OS X, still doable on windows), run volk_profile (takes a while to find the optimum SIMD implementations for your machine), and then use the following flow graph to measure resampling throughput:

It does 5 million samples per second for me. Scipy's not getting quite close to that.


The question is rather vague or inaccurate really; or I haven't understood it well.

I would start with saying "there is no such a thing as a 'best' filter (for all use)".  A filter is optimal only if it exploits a specific property of signal or noise generating better SNR over regular filter -but if you don't know about signal or noise specifically then the standard filter is also the best or optimal.

Here are some assertions that should clarify your understanding (and thereby my understanding of your question):

When we apply a filter (on any stationary signal to start with) we eliminate or attenuate part of the bandwidth and amplify the other part. So for instance if your noise is white noise (equally spread in all frequencies) and your signal is band limited to say bandwidth B, then the "optimal" filter is essentially filter that allows only the band B. If this is a baseband signal having frequencies from 0 to B, then it is a lowpass filter, if it is a modulated signal having a frequency with B width around a carrier C then the optimal filter is the band-pass filter exactly overlapping to signal.
In communication system - we transmit pulses of specific shapes. As such the pulses have a wide enough bandwidth but we know lot more specific spectrum of the signal. Here you use Matched Filter that eliminates everything that doesn't correlates to the 'known' signal.
If your noise source is very specific narrow band such as power-hum, you can use Notch Filter or Band-stop filter
One needs to consider other aspects of filter design as well. One may choose specific frequency bands for "Optimal filtering" - yet two different filters (say kaiser vs. butter-worth see this)  with same bandwidth can have different envelops of frequency spectrum - defined by order of the filter and exact impulse response. If we really want to cut down every possible frequency > B, then highest possible order filter is "Optimal" - however, if such filter creates non-uniform gain in the pass band, it may distort the signal itself while keeping noise away, so the "Optimal filter" is the one which has 'best possible' frequency response shape that produces SNR (even under the context of same bandwidth B).
Typically if your noise is not white and Gaussian, it originates from some source. You can learn more about the source of noise without signal being present to identify its property. Of course, specific noise sequence at instance t may not be of relevance at instance t + k when signal is present along with noise. But if the noise is stationary we can derive the best possible subtraction of band limited noise from a band limited signal through a theoretical formulation is called Wiener Filter. Most often in literature people are referring to wiener Filter as synonym to 'Optimal filter'.
Dealing with dynamic signals: What if either signal or noise is not stationary at all? there are host of ways dynamic noise cancellation happens - which are essentially non-stationary filters. Typically, these are adaptive algorithm which tracks noise and signal both. If you can't extract the very source of noise independently these algorithm don't quite work. Kalman filter is one such adaptive filter which works even as the system state changes dynamically. It is the most optimal filter if we know that noise is Gaussian and system state can allow signals be estimated against the actual ones. Understand a big difference between above listed filters and adaptive filters - the adaptive filters have estimates of expected signals and then measured actual signals, subtracting the two you estimate "error" i.e. noise. so through a detailed model of the system these filter already "knows" everything about Noise. If you can't estimate noise, then you can't use Kalman or other such adaptive filters. In all such stationary scenarios, theoretically nothing better can be constructed beyond the promise of "Wiener Filter".

So in essence- "What is the Optimum filter" depends on "What you know about the signal and noise"

I hope i have given right perspective with regards to your question.


How to select the specific frame with object

I'm seeing two issues here:

Frames with no activity: You can check the gradient between successive frames and only do object detection when the gradient information passes a certain threshold.
Counting the same object: The vehicle location information can be useful for this task. If the gradient change occurs in the location of a detected object in the previous frame then it is the same vehicle. If necessary, you can also use object color and object size parameters to account for varying speeds.


Your signal has the form

$$x(n)=A[\cos(2\pi fn)+i\sin(2\pi fn)]=Ae^{i2\pi fn}$$

This is a complex exponential with normalized frequency $f$ (normalized by the sampling frequency). The bins of the signal's FFT of length $N$ correspond to normalized frequencies

$$f_i=\frac{i-1}{N},\quad i=1,2,\ldots, N$$

In your case $N=8$ and $f=6/8=f_7$ (using the Matlab indexation, which starts with index $1$). So the frequency of your signal $x(n)$ lies exactly on the FFT grid. This is why there is one peak at index $7$ and all other values are $0$ (apart from numerical inaccuracies).

If $x(n)$ were real-valued, its FFT would satisfy the following symmetry condition:

$$X_k=X^*_{N-k}$$

However, for complex-valued signals, such as yours, there is in general no symmetry.


Scaling down the modal value when thresholding

Most likely this is chosen by experimental results. Seems like for the algorithm to scceed it needs a value that is lower than the modal value. But not significantly lower (this is where the 1.5 scaling factor comes in).

I suggest trying the mid-range value and its variants to see how that works. You might be able to skip the scaling procedure.
mid-range = (max value + min value)/2


First see "Naive Breakdown" in a below section. Onto PyWavelets: the algorithm was found on Github to stem from an old MATLAB implementation, but it provides no details on coding the wavelet kernel. So I set on to explore; various visuals and code annotations here. Key findings follow. -- All code
See "PyWavelets Breakdown" below. It is implementing the following - namely Eq 4:

(1) How is integrated wavelet, then differenced post-convolving, equivalent to convolving without int+diff?
The formulation is in the right, deriving from the left, equations; Instead of continuously integrating over all $\mathcal{R}$, the integral is split into segments ssumming over $k$. The input is assumed constant over $[k, k+1]$ (though shouldn't it be $)$?), as it's discrete, whereas the wavelet is kept continuous. This enables to move the input outside the integral; visually:

So how do we code this? It's clear from Eq 4 that we need a diff somewhere. pywt does np.diff(conv), which is basically np.diff(coef); why does it seem to take difference of the would-be coefficient in Naive?
In Eq 4 we have the integrated wavelet up-to k + 1, minus up-to k, multiplied with the wavelet. This is equivalent to product with the wavelet at b - 1 and b, respectively (plug in t=k+1 and t=k, compare with plugging in b=b-1 and b=b instead`).
The answer's then in in conv; assume it was already trimmed so that it's only when the wavelet is at least halfway inside the signal. Then, every point in conv is the product of the wavelet with a portion of the input, at a different tau (or rather b). So two adjacent points in conv are products with wavelets at b and b + 1:
$$
\begin{align}
  \text{conv} &= [s \cdot \psi_{\text{int}}(t - 0),\ s \cdot \psi_{\text{int}}(t - 1),\ ...] \\
  \Rightarrow  \text{DIFF}(\text{conv}) &= [s \cdot (\psi_{\text{int}}(t - 1) - \psi_{\text{int}}(t - 0)),\ ...]
\end{align}
$$
However, we seek minus between t=k+1 and t, or a greater t and lower t, which is the negative of above, thus explaining step 7 in PW Breakdown.

(2) How is resampling the integrated wavelet at increasing resolution equivalent to changing its scale?
See figure under (5) in PW Breakdown. pywt appears to "scale" via the number of samples that define the wavelet, rather than recomputing the wavelet with a dilated argument. ... though, visually, this resampling clearly does dilate the wavelet, stretching it relative to the input.
So it's indeed a form of scaling; but how does it compare to naive? Naive is recomputed more intuitively per equations, using fixed wavelet length, and passing a scaled argument. Then, which is better? The opposite question is clearer to answer; consider what happens at higher and lower scales:

Naive higher: wavelet dilates too much, exiting the sampling frame far before its zero tails.
pywt higher: wavelet is just fine; the problem's with conv. Wavelet length is fixed at 1024, so if the input is any shorter, then higher scale wavelets can never fully multiply the signal. The greater the disparity, the more the wavelet is "seen" similar to "Naive higher" by the signal; this can be seen in the question's heatmaps differing by vertical shifts.

Maximum scale of the wavelet is also unlimited; len(j) grows proportionally with scale (but there's a normalization problem, see here).

Naive lower: all's good.
pywt lower: wavelet has small resolution.
Fixed vs variable window length: higher scales must correspond to greater dilation, or being-nonzero, of the wavelet relative to the input. Both accomplish this, except Naive discontinues to zero abruptly at higher scales.

The higher scale problems in both are solved by limiting the maximum scale; if both are in their completely safe zones, the Naive implementation seems preferable due to higher resolution. However, this "safe zone" for Naive is prohibitively small, if we are to consider "safe" as not jumping to zero at tails.
To contrary, pywt's problem is more easily solved by increasing maximum wavelet length, which was suggested to rid of the zipper-like artifacts (which can be seen in the question's plots). That'd, however, worsen higher-scale performance for short inputs - but in practice it's unlikely to be shorter than 1024.

PyWavelets Breakdown:

Wavelet, prior to integration, matches exactly with the shown code blob, which is an approximation of the complete real Morlet (used by Naive) assuming $\sigma > 5$ in the Wiki.
pywt integrates real Morlet via np.cumsum(psi) * step, accounting for the differential step size
The integrated wavelet, int_psi, is reused for all scales
For each scale, the same int_psi is resampled at increasing resolution, with j as an integer varying linearly between same min and max (almost) for all scales:

The resulting wavelet for increasing scale, side-by-side with naive (for later):

The convolution result first diff'd,
negated,
multiplied by sqrt(scale),
then trimmed on both sides to either (a) only include parts where the wavelet is at least halfway "inside the signal", or (b) force its length to len(x).

Naive Breakdown: if unfamiliar with CWT/STFT, I strongly recommend parts 1-3 here.
We need (1) wavelet type; (2) wavelet length; (3) wavelet increment/overlap. We'll use real Morlet, 96 samples, and increment by 20.
win_len = 96  # "win" == window == wavelet
win_inc = 20
n_wins = (len(x) - win_len) // win_inc + 1

coef = np.zeros((n_wins, len(scales)), dtype='complex128')
kernel = morlet_kernel(win_len)

for tau in range(n_wins):
    start = tau * win_inc
    end   = start + win_len
    coef[tau, :] = _transform(x[start:end], kernel, scales)

def _transform(x, kernel, scales):
    coef = np.zeros(len(scales), dtype='complex128')
    for i, scale in enumerate(scales):
        psi = np.conj(kernel(scale))
        coef[i] = np.sum(x * psi / np.sqrt(scale))
    return coef

Visualized for scale=2:

The difference between this and naive implementation in the question is, latter uses win_inc=1, equivalent to convolving with a flipped wavelet and trimming the result to only include parts where the wavelet was fully "inside the signal". Also, here we have tau ("increment") in the outer loop, there scale, but the two are equivalent.
Lastly, shapes along timeshift dim differ; question's implementation trims so that only parts of convolution where wavelet is at least halfway "inside the input" are included; this one trims further to where the full wavelet is included.

UNSOLVED:

Why * sqrt(scale) (step 8 pywt)? Eq 4 shows 1 / sqrt(scale). I considered * scale per chain rule, except we've integrated the wavelet numerically... Makes sense to have per-scale treatment; do we get away with not integrating a per-scale wavelet because the per-scale wavelet is... the same resampled wavelet?  Resolved.


Measuring acoustic phase difference between two receivers

Keep in mind the following calculation:

$$ \frac{1130~\text{feet}/\text{s}}{1700~\text{cycle}/\text{s}} \approx 0.66~\text{feet/cycle} \approx 20.26~\text{cm}/\text{cycle} $$

Your delay distance is longer than your wavelength.  With an FFT approach, you will only get the delay as a fraction of a cycle, so you'll have to account for that.  If you can't frame your FFT on a whole number of cycles, you will need to estimate the frequency and calculate the phase shift from there.  Then, using your frequency you will need to calculate the actual delay distance.

Generally, if it is the same signal, you will get much better results using autocorrelation in the time domain.  Your next best bet is to make sure your FFT is aligned on a whole number of cycles, then you only need to calculate one bin.  If you can't do that, in my blog articles, I have frequency formulas from two or three bins and an method for determining the phase very accurately using the nearest three bins.  (See my profile page for a link to my blog).


The derivative of the Dirac delta impulse is written as $\delta'(\tau)$. This helps with notation because the mistake you made is to write $h(-\tau)=\frac{d}{d\tau}\delta(-\tau)$, which is not the case because $\delta(\tau)$ is an even (generalized) function, whereas the derivative operator $\delta'(\tau)$ is an odd (generalized) function:
$$\delta'(\tau)=-\delta'(-\tau)\tag{1}$$
Consequently, the auto-correlation function of the output of the differentiator is
$$\begin{align}R_{YY}(\tau)&=R_{XX}(\tau)\star\delta'(\tau)\star\delta'(-\tau)\\&=-R_{XX}(\tau)\star\delta'(\tau)\star\delta'(\tau)\\&=-\frac{d^2}{d\tau^2}R_{XX}(\tau)\tag{2}\end{align}$$


Determine the type of filter from the transfer function

That is an FIR filter transfer function.

Plot the frequency response to determine if it is HP, BP, LP.

(Hint: compare the magnitude at w=0 and w=pi to get a quick idea)

This seems like a homework problem, so I don't think it is appropriate to supply a complete answer.


Low-pass filter with varying cut-off frequency over time

I wouldn't call this novel but... I'm going to give an example on how to do this for a specific filter with a specific discretization.

A 1st order butterworth low-pass filter is given by:

$H(s) = \frac{Y(s)}{X(s)} =\frac{\omega_c}{s + \omega_c}$

Where $\omega_c$ is the cut-off frequency of the filter.

The backward difference (aka backward euler) approximation is given by:

$s \approx \frac{1-z^{-1}}{T_s}$

Where $T_s$ is the sampling time. By substitution, the discrete transfer function becomes:

$H(z)  = \frac{Y(z)}{X(z)} \frac{\omega_c}{\frac{1-z^{-1}}{T_s} + \omega_c} = \frac{\frac{\omega_c \cdot T_s}{1+\omega_c\cdot T_s}}{1 - \frac{1}{1+\omega_c\cdot T_s}\cdot z^{-1}}$

Let $a = \frac{\omega_c \cdot T_s}{1+\omega_c\cdot T_s}$, then:

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

Thus the difference equation is:
$y[n] = a\cdot x[n] + (1-a)\cdot y[n-1]$

Why go through the tedious process of getting an equation for $a$ when I could have just used the above formula and try to guess some values for $a$? Well, now if you want to change the cut-off frequency $\omega_c$ in real time, you simply have to recalculate $a$ each time you compute the signal. This, of course, assumes you have a constant sampling period.

Note that this filter is an IIR filter, not FIR as you requested, but the logic is basically the same: as long as you know how your coefficients relate to the cut-off frequency, then you just simply have to recalculate them every time you want a different cut-off frequency. Nothing stops you from repeating the same steps for different discretization methods, different types of filters, etc...

You may even wish to work purely in the discrete-time domain, but the relevant part is always knowing how your coefficients relate to the cut-off frequency.


If the case with two sines is clear (sine wave modulation of a sine wave carrier) then consider expanding any other shape (such as the ramp) into its Fourier components and then superimposing the result of each.
Any waveform with sharp transitions (such as a ramp) requires an infinite number of frequency components, so better to work with smoothed waveforms where the higher frequency components are essentially filtered out.
Another tip when simulating carrier modulation is to not actually simulate the carrier but use a complex baseband signal instead (what you simulate their will match what would occur at any carrier frequency and doesn't require as high of a sampling rate to simulate). This requires understanding the baseband equivalent analytic signal, but would be the recommended next step regarding understanding DSP concepts.


Is it possible to recover QAM constellation of an undersampled OFDM signal?

Well, apply the definition of OFDM to a subcarrier range of an OFDM signal: it's still OFDM, in all letters of that acronym.

So, on a superficial level, yes, this works easily¹. Looking closer:

OFDM signals aren't usually rated for some sampling rate, the systems transmitting them just work at a specific rate, which mathematically coincides with the product of subcarrier spacing and numbers of subcarriers.

The good thing about OFDM is you don't have to worry about what the spectrum looks like in between the subcarrier "frequency points"; the orthogonality guarantees you've got no signal from other subcarriers on any subcarrier frequency, and the way perfectly synchronized OFDM works is by the DFT – which essentially samples the frequency domain.

Now, have you ever wondered why in e.g. a waterfall plot or on a spectrum analyzer, the spectrum in the unoccupied carriers at the edge doesn't actually decrease as fast as the (ideliazed) 6dB / decade of sidelobe power suggest?

We often forget that yes, as you say, OFDM systems are actually discrete systems based on the DFT; often, the DACs involved are actually internally oversampling, even. This implies a "folding" imaging scenario, where sidelobe power that would fall out of the (nominal) DAC's bandwidth will be folded back – not a problem, because it just adds power from other sidelobes to sidelobes; at the transmitter, the orthogonality still stands!

Now, at the receiver, we're happy that the edge carriers are unoccupied – because the moment where we're not sampling-clock-synchronized and frequency-synchronized yet, we'll have a lot of fun with circularly shifted-in sidelobes overlaying the tentative correct subcarrier frequencies at both ends of the spectrum (assuming we critically sample).

Your subset-of-subcarrier-restricted OFDM receiver runs exactly into that hardship: No "leeway" for frequency offset and clock error to work out. Note that in this case, the frequency offset isn't "worse" to deal with than in the "normal" OFDM receiver case. It's just that you're guaranteed to have aliases from the adjacent carriers and their sidelobes also ending up in your signal – and whilst these have the same PSD as in the full OFDM signal, your processing gains on the DFT structure of the signal during synchronizations are but a fraction of what the full system has!

So, as Jason R commented:

Make sure the analog side of your receiver has adequate filtering. In theory, the sidelobes of adjacent channels don't matter, in practice, your sampling clock offset makes them hurtful. At any rate, you don't want other carriers to alias into your signal.
Make sure your receiver has plenty of dynamic range headroom. Really - OFDM has notoriously nasty PAPR, and you don't want dynamically appearing intermodulation products (hint: in many hardware SDR receivers, sampling rate and RF frequency generation are somehow coupled, so sampling rate-rationally-related intermodulation products often appear – especially bad here, because these can map "unseen" carriers to your carriers.
Use additional digital low-pass filtering to enforce unoccupied carriers at the edge of your observed band; yes, that'll take energy away from the carriers you want to see (remember, sidelobes in OFDM do carry information energy!), but it'll make your synchronization life easier.
Adapt the synchronization algorithm to your carrier subset. This might or might not be easy or even possible, considering we know nothing about how your system does synchronization.
Equalization might become a lot harder, if the original system was designed to interpolate channel state between neighboring carriers (and at bandwiths> 200 MHz, we can be pretty sure it does so to minimize the need for pilots), since you've got only a very partial view on things.

Astonishingly, the QAM aspect together with the last point might be a show-stopper here: QAM requires relatively good equalization, and since your part of the OFDM signal might give you less than an equivalent part of the pilot information in the overall system, you might simply be channel-coherency limited (either not enough pilots in frequency direction, or not enough pilots over time). So, make a channel model first, figure out how often you'll see pilots, and estimate whether things can work first, before spending too much time on this.

¹ assuming you can still frequency-synchronize. Since you said you understand OFDM very well, I assume you've got practical knowledge of OFDM synchronization methods – but in short: most do work on a restricted bandwidth, with slight modifications; you need to assume 0 subcarrier coarse frequency offset, anyways. Note that I'd still assume that if implementing an OFDM receiver, the synchronization would be the hardest part to get right. Timing synchronization is another challenge, but gets even more system-specific; if your system uses Schmidl&Cox to detect both freq and time, then you're more or less lucky and just need to adjust tolerances (and sample counts) to make this work.


How to correct for delays in a PAM-X signal with orthogonal CDMA?

Couple of observations:

you can of course also correlate with the already pulse-shaped CDMA sequence as a whole – it's a long filter, then, but that (in AWGN) will essentially give you the basis of a maximum likelihood receiver.
Commutativity of convolution: Doesn't matter whether you pulse shape filter -> CDMA sequence filter or CDMA sequence filter -> pulse shape filter. The result is the same (and the same as if you first convolved CDMA sequence and pulse shape, and the signal then with that – see 1.)

In either case, you'd only do one synchronizer.


Time domain to frequency domain conversion of audio signals to extract 1/3 octave frequencies

As per my knowledge it is desired to convert received signals from time domain to frequency domain and extract the noise levels for each 1/3 octave frequency band.

No. This is needlessly complicated. A-weighting and C-weighting can be implemented through simple FIR or IIR filtering directly in the time domain. There are public domain filter designs avialble. For example check out https://www.mathworks.com/matlabcentral/fileexchange/69-octave

  It displays noise levels in Leq sound pressure levels.

Your biggest problem will be to deal with microphone calibration. Both frequency response and overall sensitivity of Android microphones are all over the place and vary a great deal from model to model, can change between Android version and sometimes are actually actively modified by each carrier. Without some sort of calibration, you can display a number but it won't be Leq and it will mostly be meaningless.

iOS might be easier: Apple does a much better job in controlling microphone consistency in design and supply chain. For a meaningful Leq estimate you still need calibration, but it's going to be less severe and there is a much smaller number of models/iOS versions to deal with.


cancellation controller transfer function design


The error may be in calling the actual time series "expected data" (unless you expect an unpredictable random walk).

If each sample is random and uncorrelated with all previous time series measurements, except for having a distribution mean near that of the immediately preceding sample (such as in a random walk), than any predictor (NN or other) of the next sample might lag (guess near the mean of the distribution, e.g. just the previous sample).

In the cases where the mechanics of some hidden state (planetary rotational or orbital momentum, etc.) is producing the measured data, a NN might be able to be trained to deduce certain hidden oscillators (low order diffEQ solutions involving that state) or some trend line, and thus not lag.

So, the lag you find may be a measurement of how well your NN is differentiating between the two cases above.

If you don't want to see a lag, then try increasing the training error function for producing a lagged output until the NN is trained to produce a different (more random?) output than a mere input follower, given a random walk input.


If $X$ and $Y$ are zero mean independent Gaussian random variables with different variances, what is the density of $\sqrt{X^2+Y^2}$

According to your description, $x = Z \cdot \cos(W),\ y = Z \cdot \sin(W)$. Follow this answer for the derivation of joint PDF of $(Z, W)$ : Complex Gaussian Magnitude and Phase Joint PDF Derivation

You will reach the following expression after Methid of transformation:
$$f_{Z,W}(z,w) = |\mathbf J|.f_{X,Y}(z \cdot \cos(w), z \cdot \sin(w)), \ where \ \mathbf J \ is \ Jacobian$$
Since the computation of $\mathbf J$ will not change hence:
$$\mathbf J= z$$

Now, when you put the joint PDF of 2 independent Gaussian random varaibles $X, Y$ in the above expression you get the following :
$$f_{Z,W}(z,w) = |\mathbf J|\cdot f_{X,Y}(z\cdot \cos(w), z\cdot \sin(w)) $$ $$= z \cdot \frac{1}{\sqrt{2\pi \sigma_x^2}}\exp(-\frac{1}{2\sigma_x^2}(z\cdot \cos(w))^2)\cdot \frac{1}{\sqrt{2\pi \sigma_y^2}}\exp(-\frac{1}{2\sigma_y^2}(z\cdot \sin(w))^2)$$

Combine this and find the marginal PDFs of $Z$ and $W$ from the above expression. You will find that integrating from $-\pi$ to $\pi$ w.r.t. $w$ will not give you Rayleigh Distribution anymore. And, similarly, integrating from $0$ to $\infty$ w.r.t. $z$ will not give you $\frac{1}{2\pi}$, and hence not Uniform Distribution anymore.


Does Increasing Resolution Lead to Have a Poor Resolution?

Indeed, that's the Heisenberg Uncertainty Principle - you can't have both very good frequency and time resolution. You always have to sacrifice something. In case of Short Time Fourier Transform it's straightforward, but for wavelets are being 'squeezed', which is changing their frequency resolution. Figure below describes more than a thousand words:

EDIT:
Because you are still trying to understand why it is like that, then please mind that by increasing the frequency of the wavelet ($f_c$), you should also shrink it in time (decrease $f_b$). For more info please refer for example to: Paul Addison - The Illustrated Wavelet Transform Handbook. You can find there a following graphs (Section 2.12):


Help in proper notations and mathematical formulation

You observe a sequence $(s_1,\ldots, s_N)$ where random variables $S_k$ are i.i.d. and are drawn from a set (alphabet) $\mathcal{A} = \{a_1,\ldots, a_m\}$ with letter probabilities $P(a_i)=p_i$ for $1\leq i \leq m$ and $\sum_{i=1}^m p_i = 1$.

The likelihood function can be written as:

\begin{eqnarray*}
l( p_1,\ldots, p_m; s_1, \ldots, s_N ) &=& P( S_1 = s_1, \ldots, S_N = s_N ) \\
&=& \prod_{k=1}^N P(S_k = s_k) \\
&=& \prod_{k=1}^N \sum_{i=1}^m p_i \mathbf{1}\{s_k = a_i\}
\end{eqnarray*}
where $\mathbf{1}$ is the indicator function which evaluates to one when the condition is true, and zero otherwise.

In part 2 of your question, we need to change the random variable to letter counts instead of the sequence itself. Let $X_i$ be the number of times we observe letter $a_i$ for $1\leq i \leq m$ in our sequence of length $N$. Then $\{X_i\}_{i=1}^m$ has a multinomial distribution:

$$
P(X_1 = x_1,\ldots,X_m=x_m) = \frac{N!}{\prod_{i=1}^m x_i!} \prod_{i=1}^m p_i^{x_i}.
$$

(Convince yourself this is in fact the right probability mass function for letter counts; the argument involves some combinatorics.)

The log-likelihood function can be written as:

\begin{eqnarray*}
\log l(p_1,\ldots,p_m; x_1,\ldots, x_m) &=& \log P(X_1 = x_1, \ldots X_m = x_m) \\
&=& \log(N!) - \sum \log(x_i!) + \sum x_i \log p_i.
\end{eqnarray*}

The MLE of $p_i$ can be found by maximizing the log-likelihood. However, note that we have a constrained maximization problem here because we require $\sum p_i = 1$.

To convert this to an unconstrained maximization problem, let's introduce a Lagrange multiplier $\eta$. So the MLE of $p_i$ can be expressed as:

$$
\hat{p}_i = \arg \max _{p_i} \log(N!) - \sum \log(x_i!) + \sum x_i \log p_i + \eta \left( 1-\sum p_i \right).
$$

Setting the gradient of the objective function to zero gives:

$$
\frac{x_i}{p_i} - \eta = 0 \implies p_i = \frac{x_i}{\eta}.
$$

However, note that in order to meet the constraint, we must set $\eta = N$ (because $\sum x_i = N$).

This finally gives us the MLE:

$$
\hat{p}_i = \frac{x_i}{N}
$$

which can also be expressed as
$$
\hat{p}_i = \frac{ \sum_{k=1}^N \mathbf{1}\{s_k = a_i\} } {N}.
$$

This is intuitively pleasing! To get the ML estimate of $p_i$ count the number of times you see the letter $a_i$ in your original sequence and divide that by $N$.


Calculation of Reverberation Time (RT60) from the Impulse Response

Encouraged by Hilmar, I've decided to update the answer with all the steps necessary to calculate the Reverberation Time from a scratch. Presumably, it will be useful for others interested in this area. Obviously, it is the simplest approach because more advanced are definitely beyond a scope.

In the beginning, you must obtain the impulse response of a room. It can be done in various ways:

Firing a starter gun, popping a balloon, etc. - basically recording any impulsive-like signal with broad frequency content and omnidirectional characteristic. This is the simplest method of obtaining impulse response.
Sweep Sine measurement with a usage of either linear or exponential sweep. This method is my favorite as it allows you to extract many different parameters of your system at the same time.
MLS (Maximum Length Sequence) measurement with this kind of noise and using Hadamard transform to obtain the impulse response. Some comparison with sweep sine technique can be found here.

After getting the impulse response $h(t)$ of your acoustic space you are able to calculate the reverberation time and other speech related parameters such as $C_{80}$, $STI$, $D_{50}$ (Clarity, Speech Transmission Index, Definition), etc.

Very first step is to filter your impulse response in an appropriate frequency band. This is due to fact that reverberation time is a function of position in the room and frequency. Obviously, we assume that sound field is diffused and at each position in a room decay rate is the same. Thus we must still consider dependency on the frequency band. We are using 1/1 octave or 1/3 octave filters. Its desired parameters are defined by EN 61260 standard, preferably you always wish to use the filter of class 0. For most of the time, Butterworth 3rd order filter is enough (although for frequencies below 125 Hz you might experience some problems with characteristic not meeting specifications). I did the whole implementation on my own with few tweaks, but for usual applications widely used MATLAB implementation is good enough.

Next step after filtering the $h(t)$ and obtaining $h_f(t)$, is to make the response as smooth as possible before conversion to the logarithmic scale. For that purpose, Hilbert Transform is a widely used tool. The goal is to create the analytic signal:

$$h_A(t)=h(t)+j\tilde{h}(t) $$

where $\tilde{h}(t) $ is the HT of your filtered impulse response and $j$ is a complex unit. Now, knowing that your analytic signal is complex, two things are coming from its representation:

$$h_A(t)=A(t)e^{j \psi (t)} $$

where $A(t)$ is the envelope of your signal and $\psi (t)$ is an instantaneous phase. Obviously, we are interested mainly in the envelope (magnitude of the analytic signal). Below you can see overlayed filtered response and its envelope:

One clue for people speaking MATLAB'ish. The hilbert function returns the analytic signal, so no need of multiplying by $j$ and adding original $h(t)$. Just use abs(hilbert(h)) to get the envelope.

Now we have the smoothed version of our impulse response, but it must be smoothed further. It is optional and one might consider it unnecessary. It is based on Moving Average filter of appropriate length $M$. Everyone knows the moving average filter, it is basically a low-pass FIR filter with frequency response given by:

$$H(f)=\frac{\sin (\pi f M)}{M\sin (\pi f)} $$

Plots for different lengths $M$ are shown below:

And on the following figure you can observe the effect of smoothing of the previously calculated energy curve ($A(t)$ in logarithmic scale becomes $E(t)$) for two different values of $M$. For sampling frequency of $48\mathtt{kHz}$ I've used $M=5001$.

So now we have the Energy Curve (not a Power) which is simply the smoothed Hilbert envelope in logarithmic scale. Conversion to decibel scale is done by equation:

$$E(t) = 20\log_{10} \frac{A(t)}{\max A(t)}  $$

Curve $A(t)$ can be further smoothed for calculations by using the Schroeder Integration method of your envelope (also known as inversed time integration). This method gives you maximally flat decay curve and is very easy to implement.

$$L(t)=10 \log_{10} \left[ \dfrac{\int_t^{\infty}h^2(\tau )d\tau}{\int_0^{\infty}h^2(\tau )d\tau} \right] $$

Again for MATLAB people:

L(td:-1:1)=10*log10(cumsum(hA(td:-1:1).^2)/sum(hA(1:td).^2));

Mind that limit of integration (td) is equal to $\infty$. It is true when we do not have any environmental noise. Otherwise, you will see that decaying sound 'dives into noise level'. It is depicted below. The red curve is based on the correct limit, where decay crosses noise floor. If you choose your limit of integration to be too short, then as in green curve, the estimate is not long enough, so you will get under-estimation while doing linear interpolation. On the contrary, the orange curve has too large integration limit and you will get over-estimated. You always want to find the correct limit. If you didn't obtain curve in such way, then please do so.

There are many different methods for estimation of integration limit in Schroeder integral, I think that easiest to implement is the one proposed by Lundeby et al.

Regarding calculation of RT itself. It must be done by performing linear interpolation of your decay curve (or should I say Schroeder curve) with linear function: $L=A\cdot t+B$ on the correct range (which is described below). When it's done, you calculate the reverberation time from the equation. You can see that point of intercept doesn't really matter (per se), you care about the gradient of your line.

$$RT = \dfrac{-60}{A} $$

where $A$ is a slope coefficient from interpolated line (in dB/s). While you have that, you can also calculate the correlation coefficient of your linear fit (or r-value), that is telling you how good is the fit.

Finally - limits. You cannot really calculate $T_{60}$ having $60 \mathtt{dB}$ of dynamic range as you are trying to do - you always need more (unless you really know what you are doing). So you are doing linear fit on decay curve at following dynamic ranges, and the result is being extrapolated, according to ISO 3382-2:

$EDT$ (Early Decay Time): upper limit is $0 \mathtt{dB}$ and lower is $-10 \mathtt{dB}$. This parameter correlates well with perceived reverberation time. In practice, though beginning for the sake of algorithms, people are using an interval of $-1 \mathtt{dB}$ and $-10 \mathtt{dB}$ (i.e. in Norsonic analyzers).
$T_{10}$: upper limit must start at $-5 \mathtt{dB}$ to remove any fluctuations and then lower limit is taken to be $-15 \mathtt{dB}$, but it always must be at least $10 \mathtt{dB}$ above the noise floor. So in fact you need at least $25 \mathtt{dB}$ of dynamic range (or INR) to be able to calculate $T_{10}$ ($5+10+10$).
$T_{20}$: upper limit at $-5 \mathtt{dB}$, lower at $-25 \mathtt{dB}$. Minimum dynamic range needed is $35\mathtt{dB}$
$T_{30}$: upper limit of $-5 \mathtt{dB}$, lower at $-35 \mathtt{dB}$, with minimum $45 \mathtt{dB}$ of dynamic range.

Why are people using these values? Well, because you very rarely have $75 \mathtt{dB}$ of dynamic range to be able to estimate $T_{60}$. For very smooth decay these values should be equal. Some more detailed analysis of effects can be found in this publication:

  Measuring Room Impulse Responses: Impact of the Decay Range on
  Derived Room Acoustic
  Parameters

If you will not perform Schroeder integration, then most probably they will diverge a lot. For the very nasty room you can get something like:

In theory: $EDT$  = $T_{10}$ = $T_{20}$ = $T_{30}$. Although they are not necessarily going to be. For example, if your small room is coupled with large one (let's think of small chapel coupled with cathedral by a small doorway), then you will get very sharp decay for smaller one at the beginning of the curve, and very long tail in the end:

To wrap up. According to standards, you should always estimate your reverberation time $T_{N}$ starting from $-5 \mathtt{dB}$ and down to $(-5-N) \ \mathtt{dB}$. Also, ensure that you are at least $10 \mathtt{dB}$ above the noise floor. If not, then you must use another estimate. In your case, I suggest using the $T_{30}$.

Below you have example of $T_{20}$ estimation for a troublesome impulse response. $E(t)$ is the energy decay curve obtained by taking the logarithm of impulse response Hilbert envelope. $L(t)$ is the Schroeder integral. Correlation coefficient between linear fit and $E(t)$ on a given range ($-5 \mathtt{dB}$ - $-25 \mathtt{dB}$) for this case is: $cc=-0.9987$. Calculated RT is $T_{20}=1.15 \mathtt{s}$

I think that should do. If anyone will follow these steps, then he should get an accuracy of $0.1 \mathtt{s} $ against widely used Dirac. Anyway, it must be remembered that even different software tend to give different results (i.e. EASERA and Dirac), so you should be totally fine. For more complicated impulse responses with artifacts at the end of impulse response one might get worse performance, but anyway, such measurements are not reliable and reflect wrongly conducted an experiment.


What are the optimal features for text line and table row discrimination in a handwritten document image?

You can detect the table from a normal text paragraph using the Hough Transform by which you can differentiate the text from the image. You can also use the bounding box method to detect the table and differentiate from the normal text rows.


Yes it is. However, unlike real AWGN which comes mostly from thermal agitation in electronics, matlab's noise is not purely random, this might be the cause of you seeing a deterministic pattern. But this pseudo-random noise is usually random enough to approximate real life noise

If you want to assess the noise "randomness", plot an intensity histogram of In-I, with In being the noised image and I the original image (you are subtracting the deterministic part and keeping the random part).


You just need to split the sum into the positive and the negative indices (without forgetting the index $k=0$):

$$x(t)=a_0+\sum_{k=1}^{\infty}a_ke^{jk\omega_0 t}+\sum_{k=-\infty}^{-1}a_ke^{jk\omega_0 t}\tag{1}$$

By simply changing the sign of the index $k$ in the second sum in $(1)$ you get

$$\begin{align}x(t)&=a_0+\sum_{k=1}^{\infty}a_ke^{jk\omega_0 t}+\sum_{k=1}^{\infty}a_{-k}e^{-jk\omega_0 t}\\&=
a_0+\sum_{k=1}^{\infty}\left[a_ke^{jk\omega_0 t}+a_{-k}e^{-jk\omega_0 t}\right]\end{align}\tag{2}$$


How to find phase values on precise pixel locations of a given Image?

You can't find the precise phase value of a given pixel because of the uncertainty principle. You can only find the phase value for the pixel's neighborhood using short-time Fourier tranform.

If you want to know the phase for pixel (100,100), ou have to take a neighborhood, for instance (50:150,50:150) and calculate the 2D Fourier transform on this window.


What are the various techniques for detecting Walls in the Building(Architectural) floor plan Images?

Try using Template matching by taking a small cutout of the wall whose length is long enough so that it is not confused with other elements. Keep the threshold value high(around 0.9 or higher) and then check the bounding boxes predictions.


How to find correlation/cross-correlation of two signals in real time?

Correlation is a multiply and accumulate process. So to do what you describe, simply multiply the two signals and accumulate (integrate) the resulting output. The integration time is set based on your knowledge of how long the waveforms could be correlated; a longer integration time will give you more processing gain against background noise but will also effect how rapidly you can observe changes in the correlation. If the integration time is longer than the interval over which the correlation changes the resulting correlation will be averaged out (therefore not observable).

Here is a specific approach that would minimize processing required while giving you a new result after every incoming sample, when we can assume synchronous inputs:

Multiply the two inputs, and place the result in a "sliding FIFO buffer" (FIFO = First In First Out). The buffer is of length N where N is the number of samples in our "integration time". This is only for purpose of maintaining what the resulting multiplication was N samples ago, as there is no need to do operations on every element in the buffer on each update or to keep a buffer of the input samples.

For each input, multiply the two and add the resulting product p to the buffer, which causes the last sample in the buffer to "fall off".  Add the result to an accumulator register r while subtracting the product from N samples ago (now in the last position in the buffer):

$r(n) = p(n) - p(n-N) + r(n-1)$

Where:

$r(n)$ is the current value to store in the accumulator

$p(n)$ is the product of incoming samples $x(n)$ and $y(n)$

$r(n-1)$ is the previous value that was in the accumulator before this update

$p(n-N)$ is the value of the product from N samples ago (get this from the FIFO buffer).

Thus at any given time r(n) represents the (unnormalized) correlation of x and y over N samples.

Just to note, and not important if you don't understand this part: This is a CIC (Cascade Integrator Comb) implementation of a moving average filter. A block diagram showing this process is shown below:

The above will give us a result that is not normalized. If a true normalized cross correlation was needed, we would need to scale the result by the number of samples as well as the standard deviation of x and y. Note that same approach can be used to determine the variance for both x and y: Use the same block diagram with x as both inputs to determine the variance of x (since the variance is the sum of the squares divided by N), and similarly y as both inputs to determine the variance of y. In most cases we would want to avoid this if we can be content with a result that is proportional to the correlation given the increased processing requiried due to the need to divide and do square roots (and three of the above block diagrams). If you have to go down this path, look into the Cordic Rotator for computing the square root if you can spare doing multiple operations between each incoming sample, and use an N that is a power of 2 so that the divide by N is a simple bit shift.


How to reduce frequency resolution for high sampling rate and lot of samples?

A simple way is to split the data into smaller N sample segments (pick N to suit your desired frequency resolution, then average the power or amplitude of all the FFT spectra:

$$
  \hat X = \frac{1}{N} \sum^{N-1}_{n=0} X_n  X_n^*
$$

Where $X_n, n \in \left\{0:N-1\right\}$ is your series of shorter FFT results and $\hat X$ is the averaged power spectra.

If I remember correctly, for noise signals, the variance is reduced by $\frac{1}{N}$.


How do I implement LIBSVM for multiclass classification using java?

libsvm doesn't work with only integer values -- in fact, you're supposed to scale your features to lie within [-1,1] or [0,1] before applying SVM's typically. Read their guide. They also have some tricks for dealing with categorical features as well.

If you want to get finitely many values which the feature can take, use quantization.


Downsampling fft results to a graphics display

Ok, here is the answer...hopefully this will help someone else.  Read this thesis:

http://skemman.is/stream/get/1946/15343/37285/3/SS_MSthesis.pdf

This gives a very easy and comprehensive overview of how to do just what I needed. I do believe this particular case is different than the usual downsampling. It fits the situation perfectly.
The thesis describes several methods of downsampling for display purposes. For example the practice of just dropping every nth sample, performing an average over the bucket of samples corresponding to the downsample ratio and more advanced methods.  The one I will be looking at is called Largest Triangle Three Bucket.
Basically you divide your original samples into buckets consisting of samples equal to your downsampling factor.  For example if you have a sample set of 10,000 points and you want to only display 2000, you will divide your original samples into 2000 buckets containing 5 sample each.
You then rank every point in the bucket by calculating the area of a triangle it forms
with the selected point in the last bucket and the average point in the next
bucket.
This seems to work really well in preserving the visual representation of data.


Argument of Dirac delta function is quadratic with complex roots

I shall provide more details if I am correct. IMO if there is not root on the domain of integration, and here I suppose that $t\in \mathbb{R}$, then the argument never vanishes. Then, in an engineer fashion, one could say (but $t\mapsto \delta(t)$ is not a function):

$$\delta(t^2+t+1) \equiv 0\,.$$

In more precise words, I would consider that for any suitable $f(t)$, $t\in \mathbb{R}$:

$$\int_{-\infty}^{\infty}f(t)\delta(t^2+t+1)dt = 0\,.$$

Related links:

Dirac $δ(g(x))$ with complex roots $x_i $
Regarding $x^2−a^2$ inside the argument of Dirac delta


Is this calculation of inverse z-transform proper

HINT:

Write $H(z)$ as

$$H(z)=G(z^2)\tag{1}$$

with $$G(z)=\frac{z^{-1}}{1-0.5z^{-1}}\tag{2}$$

Determine the inverse transform $g[n]$ of $G(z)$. Figure out what replacing $z$ by $z^2$ means in the time domain. From this, the inverse transform of $H(z)$ is easily obtained from $g[n]$.

EDIT: Your partial fraction expansion lacks the factor $z^{-2}$, and the major mistake is that you transform both terms to the same sequence, even though the poles have opposite signs. If you do it right, the two terms of the sequence only cancel each other for odd $n$, and you can combine them for even $n$. The method shown as a hint in this answer is probably more straightforward because it avoids partial fraction expansion.


Or $2\pi$, or $3\pi$, or any integer multiple of $\pi$.  Any odd multiple corresponds to -1 + 0i and any even multiple corresponds to 1 + 0i, aka -1 and 1.

"Phase of a real number" is a little bit of a misleading label.  What is required here is an understanding of the complex plane and what "phase" means in terms of a DFT bin value.

Your question is equivalent to "For what values of arg(z) is z a real number?"

If that is meaningless to you, I suggest you start by reading two blog articles of mine:

The Exponential Nature of the Complex Unit Circle

And the newest:

Angle Addition Formulas from Euler's Formula

There are of course many other searches.  Your terms should be "complex plane real values" for a start.

This is essential foundation material for a lot of DSP concepts.


Reconstruction of Contiuous - Time Signals

Because sampled signals Spectrum will have copies of original spectrum at multiples of $f=F_s$. Even if aliasing happened due to choice of $F_s$, the spectrum of sampled signal $x(n.T_s)$ will be periodic in frequency always with period $F_s$. And you can take any one such $k^{th}$ copy or period and integrate from $\frac{k-1}{2}F_s$ to $\frac{k+1}{2}F_s$. Now if you vary $k$ from $-\infty$ to $\infty$, then it is equivalent to integration on LHS.

Key: Sampling at $F_s$ will create a spectrum which is periodic in frequency with period $F_s$.


FFT returns large low-frequency power - probably because signal is interpreted as containing cycles that do not fit. How is this called?

Your reasoning sounds correct. It is impossible to see from your example what frequency that low frequency is.

If the frequency of the phenomenon is 0 Hz, then it is simply a direct current (DC) bias that affects the first bin of the discrete Fourier transform (DFT). If none of your signal of interest is in that bin then you can simply ignore it in your analyses or set it to zero. If your signal never goes to negative values, then you have DC bias (although there could be other low frequencies).


First of all - I myself am not a pro in Control theory, but a mathematician - so I write what I think might be what your prof means.

Part where i am quite certain:

Not to be polynomial means, that there exists NO polynomial which describes the transfer function. E.g. let $H$ be your transfer function, then there exists no $n\in\mathbb{N}$ and $p\in P^n[X]$ such that $$H(s) = \sum_{i = 0}^{n}p_is^i$$
holds.
Examples for functions which are not polynomial are: exponential functions (e.g $e^s$) or trigonometric functions (e.g. $sin(s)$), but there are many more.

The next part is where I am not 100% sure, because it is not my main subject.

I assume you mean the inverse Z-Transform when saying anti transform. I do not know too much about the z-transform, but i think it is the discrete-time equivalent to the laplace-transform and therefore:

Let $H$ be a transfer function. Then $H$ is the Laplace-transform of an impulse response. Thus the inverse Laplace transform (i.e. $h$) exists. (The stability however is something different.)

But if you want to check for a given "transfer function" (an arbitrary function in the "s"-Domain) if there exists an impulse response corresponding to this "transfer function", then I think you need the inverse Laplace transform to converge, i.e.
$$ \frac{1}{2\pi i} \lim_{y\rightarrow\infty}\int_{\gamma-iy}^{\gamma+iy}e^{st}H(s)\,ds$$ to exist for some $\gamma\in\mathbb{R}$.
The integral then is your impulse answer. There exist non-polynomial functions for which this integral exists. Thus, there are non-polynomial transfer functions which correspond to a system. Many many authors in literature however assume only rational valued transfer functions... i really do not know why.

If this answer is NOT correct (because for the discrete-time case there is something mayor different), or does not help, please say so and apology. As I said, this is not my main subject.


Your image with increased contrast:

Figure 1. Denis Kulagin's figure with gain applied.

That is kind of color banding due to quantization. Quantization errors can be made worse by repeated decoding and encoding, because quantization levels and things that are quantized may not match between the encodings. In video encoding, it is not directly the pixels that are quantized, but things like discrete cosine transform coefficients for 8x8 or so blocks of chroma and luma values of pixels. JPEG images compressed with a low quality setting exhibit color banding for the same reasons:

Figure 2. Top: 60 KB PNG file, Bottom: 2 KB JPEG file

Also color space conversions may introduce quantization errors and banding. For example sRGB and Rec. 709 transfer characteristics differ at worst by about a factor of two at low intensities.


Line tracing folowing a path of an almost connected components

Here is a possible solution:

Use FFT2 to find angle of text and approximate spaces between lines
Rotate the image to be at 0 angle
Sum the image column-wise
Find threshold that splits lines of text from non text lines
Show lines on rotated image.

It can be improved by finding a more complex curve (rather than line).

Here are the results:

Original

After rotation

Sum of rotated image column-wise

Final results

And the code is:

%% 1) Read image
im = imread('http://i.stack.imgur.com/zbHTp.jpg');
im = im(:,:,2);
im = double(im < 100);
im = padarray(im,[40 0]);
figure();imshow(im);

%% 2) Find approximate line spacing
F = fftshift(fft2(im));
% figure;imagesc(log(abs(fftshift(fft2(im)))));
[~,iMax] = max(F(:));
[r1,c1] = ind2sub(size(F),iMax);
F(iMax) = 0;
[~,iMax2] = max(F(:));
[r2,c2] = ind2sub(size(F),iMax2);
lineSpaceY = size(im,1)/ abs(r2-r1);
lineSpaceX = size(im,2)/ abs(c2-c1);
angleDeg = abs(90-rad2deg( atan2(lineSpaceX,lineSpaceY)));
imr = imrotate(im,angleDeg);
figure;imshow(imr);

%% 3) Take a look at 1D signal
S = mean(imr,2);
isText = S > max(S)/2;
isText = imclose(isText,strel('line',lineSpaceY/10,90));
% figure(); plot(isText,'o-');
rp = regionprops(isText,'Centroid');
centerY = cat(1,rp.Centroid);
centerY = centerY(:,2);

%% 4) Show final results
figure();imshow(imr);
for i=1:numel(centerY);
    hold on;
    plot( [0,size(im,2)],[centerY(i) centerY(i)],'g','LineWidth',2);
end


Inverse Discrete Time Fourier Transform of a given Function using the definition

Yes, you can do that. Because we substitute $e^{j\omega}$ for $z$ to get the frequency response of a discrete system, you can go back by substituting $z$ for $e^{j\omega}$. This yields

$$H(z) = \frac{1}{(2-z)(2-z^{-1})} = \frac{1}{5-2z^{-1}-2z} = \frac{z^{-1}}{-2+5z^{-1}-2z^{-2}}$$

Therefore the system you're trying to find is

$$-2y[n]+5y[n-1]-2y[n-2] = -x[n]$$

$$y[n] = \frac{-x[n]-5y[n-1]+2y[n-2]}{-2} = \frac{1}{2}x[n] + \frac{5}{2}y[n-1] - y[n-2]$$


Resonant peak frequencies for phaser effect

I refer to the following structure of a phaser (from this article):

Let's call the depth factor $a$, and the feedback factor $f$ ($|f|<1$).
I assume that there are $2N$ first order allpass filters, each with transfer function

$$A_0(z)=\frac{z^{-1}-b}{1-bz^{-1}},\quad |b|<1\tag{1}$$

The allpass frequency response is obtained by evaluating (1) on the unit circle $z=e^{j\omega}$:

$$A_0(e^{j\omega})=\frac{e^{-j\omega}-b}{1-be^{-j\omega}}=
\frac{e^{-j\omega}-2b+b^2e^{j\omega}}{|1-be^{-j\omega}|^2}=\\
=\frac{(1+b^2)\cos(\omega)-2b+j(b^2-1)\sin(\omega)}{|1-be^{-j\omega}|^2}\tag{2}$$

And since $A_0(z)$ is an allpass function

$$A_0(e^{j\omega})=e^{j\phi(\omega)}\tag{3}$$

From (2) and (3) the phase $\phi(\omega)$ of the allpass satisfies

$$\tan(\phi(\omega))=\frac{(b^2-1)\sin(\omega)}{(1+b^2)\cos(\omega)-2b}\tag{4}$$

If you cascade $2N$ of these allpass filters, the frequency response of this cascade is

$$[A_0(e^{j\omega})]^{2N}=e^{j2N\phi(\omega)}$$

Let us first consider a system without feedback, i.e. $f=0$. The total frequency response of such a system is

$$H(e^{j\omega})=1+ae^{j2N\phi(\omega)}=
1+a\cos(2N\phi(\omega))+ja\sin(2N\phi(\omega))$$

and its squared magnitude is

$$|H(e^{j\omega})|^2=1+2a\cos(2N\phi(\omega))+a^2\tag{5}$$

From (5) it is obvious that the locations of the maxima of $|H(e^{j\omega})|^2$ are independent of the depth paramter $a$, and that they only depend on $N$ and on the allpass phase response $\phi(\omega)$ (i.e. on the paramter $b$). The maxima occur where the argument of the cosine in (5) equals $2\pi k$, $k=0,1,\ldots,N$, which is equivalent to

$$\phi(\omega_k)=\frac{k\pi}{N},\quad k=0,1,\ldots,N$$

where $\omega_k$ are the peak frequencies. With (4) this condition becomes

$$\frac{(b^2-1)\sin(\omega_k)}{(1+b^2)\cos(\omega_k)-2b}=\tan\left(\frac{k\pi}{N}\right),\quad k=0,1,\ldots,N\tag{6}$$

Equation (6) can be used to determine the peak frequencies $\omega_k$. Note that this needs to be done numerically since this equation cannot be solved for $\omega_k$ analytically. There are just some special indices $k$ for which an explicit solution is possible:

$$\begin{array}{rcl}\omega_0&=&0\\
\omega_{N/2}&=&\arccos\left(\frac{2b}{1+b^2}\right)\quad\textrm{if } N \textrm{ is even}\\
\omega_N&=&\pi
\end{array}$$

The following figure is an example for $N=4$, $a=0.7$, and $b=0.6$ which shows $|H(e^{j\omega})|$
and the locations of the maxima (shown by 'x') which I've computed using equation (4).

If there is feedback ($f\neq 0$) then the frequency response looks more complicated:

$$H(e^{j\omega})=\frac{1+(a-f)e^{j2N\phi(\omega)}}{1-fe^{j2N\phi(\omega)}}$$

and its squared magnitude is

$$|H(e^{j\omega})|^2=\frac{1+2(a-f)\cos(2N\phi(\omega))+(a-f)^2}{1-2f\cos(2N\phi(\omega))+f^2}$$

Luckily it turns out that feedback doesn't change the location of the peaks, it just makes them 'peakier'. In the following figure the system the frequency response shown above is compared to the response with feedback ($f=0.5$):


Confusion regarding lacunarity of koch curves with same fractal dimension but different orders

I think this is a conceptual issue. Your images are not fractals, they are approximations to fractals. Lacunarity is a property of the fractal curve, not of its approximation. The Koch curve has a specific lacunarity.

Your program computes a measure that is inspired by lacunarity, from a given approximation to a fractal. This measure determines what fraction of the image can be considered a hole at a given scale.

Your normalized lacunarity I hadn’t seen before. But it normalizes the curve by the value at r=1, the smallest possible scale in the image. Obviously at this scale the difference between your two images is the largest, since at a larger scale they should be identical (the difference is in the fine details). By normalizing by this value, it is expected that the ratio elsewhere, where the differences are smaller, will reverse.


A couple of confusions here.

Features refer to some form of lower dimensional descriptors that explain a (potentially local) region of interest. They are useful in converting the appearance into certain signatures that are easier to handle / more robust than naively using pixels.

What you are trying to ask is probably keypoint detection. This term is used to refer to certain points of interest in the image, that are interesting or worthy to extract features at. Indeed, this is a more general term for corner detection. Most of the time (especially nowadays in the context of deep learning), keypoints need not be some geometrically salient corners but rather points that are critical in describing the image/signal as a whole. The state of the art methods learn this notion of saliency and critical-ity.

In fact, the term feature detection should probably refer to the combination of two steps: first detecting the keypoints and subsequently extracting descriptors on those locations.
Hope this makes sense.


Use of Long Training Field (LTF) of 802.11 standard

I'm not to familiar with the HT-LTF, but OFDM channel est. is done as follows:

rx'd long training symbol ($LTF_{rx}$) is the tx'd LTF ($LTF_{tx}$) convolved with the channel:

$$LTF_{rx}(n) = \sum_{m=-\infty}^{\infty}h(m)LTF_{tx}(n-m)$$

translate this to frequency domain:

$$LTF_{rx}(w) = H(w) LTF_{tx}(w)$$

and the channel inverse is then:

$$H(w)^{-1} = \frac{LTF_{tx}(w)}{LTF_{rx}(w)}$$


One first suggestion for you is to implement DFT filter banks in some language (Octave, Python, Matlab) and play with them to get some insight. Then,
let me provide a information source, and a few hints.  In  The DFT Filter Bank, you can find:

  In this section, we will show how the DFT can be computed exactly from
  a bank of $ N$ FIR bandpass filters, where each bandpass filter is
  implemented as a demodulator followed by a lowpass filter

The above website is really worth reading. You can see for the plot that you have two responses, one being obtained via a shift in frequency for the other. Hence you know that, from one of the Fourier modulation properties, that a shift in frequency corresponds to a modulation in the time domain. Let us now look at the low-pass filter only (in blue). A running low-pass filter has a frequency response related to the sinc function. More precisely, it is called a periodic or an aliased sinc, with frequency response:

\begin{eqnarray*}
H(e^{j\omega}) &=& \frac{1-e^{-j\omega N}}{1-e^{-j\omega }}
= \frac{e^{-j\omega N/2}}{e^{-j\omega /2}}
\frac{\sin(\omega N/2)}{\sin(\omega /2)}\\
&= &
Ne^{-j\omega(N-1)/2} \hbox{asinc}_N(\omega ).
\end{eqnarray*}
where $N$ is the length of the running sum. It can be inferred  from the number of ripples, corresponding to a number of zeroes of the aliased sinc.
The number of filters is potentially trickier, because of oversampled filter banks, but this may be outside the scope of your lecture. Look at how many shifted (modulated) versions of the low-pass can cover evenly the spectrum.

Then again, do not hesitate to verify via simulations.


As far as I can tell from the graph, the variance of the signal goes up substantially under "oscillation" conditions. So, monitor the variance over a rolling window. High variance indicates oscillation. To choose the window width, consider:

if the window is too short, the computed variance will be too noisy
if the window is too long, the monitor will be slow to respond

Also note that, given the information you've provided, there is no way to tell the difference between the step function signal itself and the beginning of an oscillation condition. You should therefore impose some auxiliary condition on the signal to avoid false positive detection. For a square wave signal, I'd suggest detecting extremely large time derivative. If the derivative exceeds some critical value, you probably have a square wave step, not an oscillation. So the criterion is:

oscillation if variance exceeds $V_{min}$ and max derivative is smaller than $D_{max}$.


How to remove the modulation before doing frequency offset estimation?

For lower-cardinality PSKs like the QPSK, the "classical" way is to
take the signal to the $M$th power, $M$ being the number of constellation points.
From the shape of an $M$-PSK modulation, it's clear that the constellation points are $e^{j2\pi \frac mM},\, m =0,1,\ldots,M-1$, and putting that to $M$th power yields
$${\left(e^{j2\pi \frac mM}\right)}^M=e^{j2\pi m}=e^0=1,$$
due to the $2\pi$-periodicity of the complex exponential.
There you go, undergrad textbook-level data removal :)
Since your frequency offset is but a multiplied $e^{j2\pi f_\text{offset}t}$, you get a constant tone at $M$ times the offset frequency after.
Small problem: you don't have only signal, you do have noise in your received signal, which gets put to the $M$th power, too, and also intermodulated with the signal (you're converting nice additive noise to multiplicative noise). That typically drastically reduces your SNR – you don't get the most efficient frequency estimator that way. I'd advise against this whole approach!
The typical way to counter frequency offset in benign-sized modulation is to do a second-order PLL with loop filter tuned to the symbol rates you're expecting.
Add in that you'll of course have Gray mapping and forward error correcting codes, and that you can thus use your corrected decisions to refine your frequency error estimate, and you'll quickly see how undesirable "deleting" the info through nonlinear means (which are the only means to do that) is. For BPSK, the squaring approach is commonly done (because you'd want to do something similar to recover the timing, not only the frequency error, anyways), as it's low-complexity.
A second-order PLL works relatively well (even abused for QAM) until your constellation becomes large, where the probability that your PLL either "swallows" a legit symbol transition or misinterprets noise become problematic.
As everywhere, the synchronization of larger PSKs quickly becomes very system-specific and usually reflects the assumptions made on the channel. For example, if you have a satellite channel with substantial Doppler, and Doppler rate, but with low phase pertubations other than that, periodically sending a synchronization-aiding words from a smaller PSK might be a very good tradeoff in terms of complexity, robustness and data rate loss.


Derivation of Nyquist Frequency and Sampling Theorem

It really boils down to aliasing.  In continuous-time, if you have any two signals $x_1(t) = \sin(2 \pi F_1 t)$ and $x_2(t) = \sin(2 \pi F_2 t)$, then as long as $F_1$ and $F_2$ are distinct, the signals are, too.

But consider sampling at some time interval $T_s$, so that the sampled signals are $x_1(k) = \sin(2 \pi F_1 T_s k)$ and $x_2(k) = \sin(2 \pi F_2 T_s k)$.  If you have a condition where $2 \pi F_1 = 2 \pi F_2 \pm 2\pi n$ for any integer $n$, then the signals are not clearly not distinguishable, by the simple trigonometric identity $\sin \theta = \sin \theta + 2\pi$.  Do the math, and you find that the frequencies cannot be separated by more than the sampling rate to be distinguishable.

It's worse, though, because you also have the trigonometric identity $\sin \theta = -\sin \theta + \pi$.  When you consider that the only thing that allows you to distinguish one signal from another is the frequency, that means that if $F_2$ is a mere $\frac{1}{2 T_s}$ away from $F_1$, then $x_1(t)$ cannot be distinguished from $x_2(t)$.

This phenomenon of two distinct signals appearing as one after sampling is called aliasing.  It's pretty much the basis of the Shannon-Nyquist sampling theorem, and the Nyquist rate.  Because if you're sampling at $\frac{1}{T_s}$ and trying to distinguish signals further apart than $\frac{1}{2T_s}$, you just can't.

Hence, the Nyquist rate.


You're doing interpolation "wrong", at least in the sense I assume you want it to happen.

So what you'd normally do is just take every input sample, and insert 15 zeros in between each of them. Spectrally, the resulting signal would be your input signal, plus 14 images of it, spaced 1/15 of the output bandwidth.

Now, your approach of sample and hold interpolation is the same signal, i.e. the zero padded one, unshifted, plus shifted by one sample, plus shifted by two..., plus shifted by 14 samples. So let $x$ be the input $y$ be the zero-padded and $w$ be your repetition-interpolated signal:

$$\begin{align*}
x[n]&:=\text{input}\\
y[n] &= \begin{cases}
x[\frac{n}{15}]& \text{for}\, n\, \mathrm{mod}\, 15 = 0\\
0& \text{else}
\end{cases}\\
w[n] &= x\left[\left\lfloor \frac{n}{15}\right\rfloor\right]\\
&= x[n] + x[n-1] + \dots + x[n-14]\\
&= \sum_{l=0}^{14}x[n-l]
\end{align*}$$

Hence, the spectrum of $w$ is just that of $y$, superimposed with 14 phase-rotated versions of itself.

$$\begin{align*}
W[z]
&= \sum_{l=0}^{14}\mathcal{Z}\left\{x[n-l]\right\}
\end{align*}$$

In fact, you've built a signal whose 15 polyphase decomposition components are identical!

Now, that will lead to interesting effects. In contrast, if you did the usual interpolation, i.e. took my $y[n]$ and low-pass filtering it to its original Nyquist bandwidth (i.e. $\frac{1}{15}$ of the output rate) it, none of the  polyphase components of the resulting signal wouldn't be identical to the previous one -- they'd incorporate a phase shift that ensures continuity.


is any filter whose impulse response has a discontinuity in the first derivative fundamentally high-pass?

This actually isn't the case. Consider the continuous-time triangle function:

$$
\text{tri}(t) = \begin{cases}
1 - |t|, && |t| < 1 \\
0, &&\text{otherwise}
\end{cases}
$$

This function has the Fourier transform:

$$
\text{Tri}(f) = \text{sinc}^2(f)
$$

This follows readily from the fact that the triangle function can be constructed by convolving a rectangle function, which has a Fourier transform of $\text{sinc}(f)$, with itself. The $\text{sinc}$ function decays to zero as its argument approaches infinity (and its square decays even faster yet), so the triangle window does not have a highpass characteristic.


Using MFCCs for acoustic machine failure prediction

The biological foundation for MFCCs are pretty weak. The mel-scale, which the mel-spectrogram and thus MFCC uses, is a rough approximation of how frequency differences are perceived along the frequency axis. The log transform of amplitudes is better than linear, but not really a good model of the ears loudness perception - which is non-uniform both with amplitude and frequency, as well as state-dependent (temporal masking).
The choice of hop-lengh is also maybe roughly close to what the human ear has, but vastly simplified - the ear temporal resolution is not uniform either. And no binaural representation either, a critical part of the human auditory system.
The primary reason for using MFCCs is its numerical convenience. The log mel-spectrogram (when computed with appropriate parameters) preserves a large amount of the information in an audio signal - but reduces the feature space considerably and introduces useful structure (time-frequency). Using a 32 ms hop, 512 samples at 16kHz = 32 frames per second, and 64 mel bands one can reduce one second of audio from 16000 dimensions to 64x32=2048. Then when computing the DCT-II to get the MFCC, most of the energy of this ends up in the lower coefficients. This lets us drop the higher coefficients without loosing much information. Typical would be 20 (or 13) coefficients, leading to 20x32=640 dimensions. The information in the different coefficients also tends to be much less correlated than the input spectrogram, an advantage for some models - such as the commonly used GMM-HMM (Gaussian Mixture Model, Hidden Markov Model).
Bottom line is that the feature space with 640 MFCCs per second, instead of 16000 samples per second, is much easier to learn. This dimensional reduction and feature transformations of MFCC is very useful, even if the input is not speech. So one can make OK performing systems for many non-speech tasks using MFCC.
That said with the advent of deep learning and larger amounts of data and compute - MFCC based models are being outperformed by log mel-spectrograms processed by 2d Convolutional Neural Networks, as well as raw waveforms processed by 1d CNNs.
Note also that in the majority of cases, we are still struggling to achieve human level performance. And with supervised learning, our labels are typically provided by humans - a performance roof which the model cannot easily surpass. In this scenario it is not that problematic that our audio reception is limited by loose inspirations from the human auditory system.
There are notable exceptions though, such as ultrasound (for bat detection etc) and infrasound (for condition monitoring) - there it is critical to go beyond the frequency range that human hearing has.


Can I use independent component analysis if I have two discrete populations of signal sources mixed together?

Conceptually, for ICA to work, it requires N observations of the SAME phenomenon and it will attempt to decompose these signals (the mixtures) to < N components.

From your description each detector is counting a different energy particle. That is, the "Mixtures" are actually already three different signals. Consider the example where you use the recordings of 4 different microphones to try and distinguish between 2 different (mutually independent) sound sources. In this case each microphone is recording a (different) sum of BOTH sound sources. In your case, because of the phenomenon you are observing, each "microphone" (detector) is NOT recording a (different) sum of BOTH sources, but rather it records each source separately (most of the times).

As you describe, sometimes, the SAME event might be recorded by more than one detectors (spill-over). In this case, if you insist with ICA then it would attempt to decompose THAT sum.

Therefore, perhaps you should consider something simpler, like obtaining the cross correlation between the events recorded by each one of the three different detectors as an indicator of the amount of spill-over.

But in this case, instead of examining 1000 "triplets" of observations, you create three different time series. One time series per detector describing the events (in time) of each detector.

If your detectors and aparatus are fast enough, this signal will look like a series of exponentially decaying pulses. If not, your signal will look like a train of pulses convolved with the impulse response of the system (a smoothed version). In either case, if the same event has been detected by more than one detectors then the output of these detectors will simultaneously rise in level (and increase the cross - correlation or some other form of correlation metric). It will be easier to detect that rise in the first case (fast system).

Therefore, the problem now becomes one where you have to decide which detector to assign that event to for it to be counted correctly. There is a number of ways that you can do this (not using ICA) but you will have to know the spectral response of your detector (Pulse amplitude VS energy or Voltage output VS energy). You would then be able to create a weighted sum of the counted particles by looking up the "weight" for each pulse on your detector's response. Or you could apply a "window" to the pulses. That is, consider a sharp threshold around the peak response of the detector and only count events that generate such (usually high) pulses. Considering a short region around the peak performance of the detector will reject the (lower) pulses created by spill-over.

I hope this helps. I hope to have understood your configuration correctly.


Sampling: How should I sketch a signal spectrum?

Figure 7.4 just concerns itself with an example signal, which happens to have triangular spectrum.

your $x(t)=a\cos(2\pi f t)$ signal is just that: a single cosine tone, and thus has the "shape" of two dirac impulses with integral area $\frac a2$ each.


Posted for anyone who may find this useful...

I created a picture that shows DFT frequency bin spacing for odd and even cases of N where N is the number of samples.  FFTs usually operate on an even number of samples (the algorithm works by repeatedly breaking the problem into halves), so only the even case applies.  The DC component (0*fs) is always part of the computation .  The bin at the nyquist frequency (fs/2) only shows up for the even case.  The frequencies above the nyquist frequency reflect back on the negative side of the frequency axis in reverse order.  For real signals the reflected FFT values are complex conjugates of their positive counterparts.  There are two specific cases depicted on the complex plane, N=5 for the odd case and N=8 for the even case.


I seriously doubt that there is a simple way to accomplish this.  Most of the work that I have seen tends to go in the direction of improving the synthesis, or using better (and more) phoneme recordings to get better quality.  If there were a simple way, I think someone would have found it by now rather than putting so much more effort into better synthesizers.

That said, you could try using a bandpass filter on the output.  The sample you referred to sounds better (at least is doesn't buzz) if you filter out everything above 800Hz, and everything below 100Hz.  That does kind of make it sound like the person speaking has a blanket stuffed in to her mouth, though.

If you would like, I could put together some python code (though you will need numpy/scipy for it) or you could use audacity or sox to see if the results are even acceptable before implementig code.


Is it possible to get the envelope of a signal via FIR filter?

Getting the envelope of a signal requires more than a simple band-pass filter.

However, you could rectify the signal and bandpass-filter the rectified signal to get the enveloppe. That's how I used to do it with ultrasound signals.

https://en.wikipedia.org/wiki/Envelope_detector

You could also use more elaborate techniques, like the Hilbert transform (or Hilbert filter) to get a complex signal that would allow use to use the modulus (or magnitude) of the complex signal to get the envelope.


A Simple Algorithm to Filter / Smooth / Denoise a Noisy Staircase Graph

Assuming you meant to produce something similar to the green line:
What about
$$\text{output}[n] = \max\{\text{input}[n-k], \text{input}[n-k+1], \ldots ,\text{input}[n]\}$$
i.e. you just find the maximum along a sliding window over the last $k$ input values?


How we calculate Precision-Recall Curve?

Obviously, you can't plot a curve from a single pair of values, so you have to find a parameter to vary that changes your problem in a way that makes sense to what you're trying to do.

For example, suppose your binary classification problem has a choosable threshold.  If you vary the threshold up and down, then it will change the precision and recall.

Then you might get a plot like the one below.  R code below that generated it.

R Code

#26568

Ndata <- 1000

c1 <- runif(Ndata,1.5,3)
c2 <- runif(Ndata,1,2.5)

data <- c(c1, c2)

thresholds <- seq(1,2.99,0.01)

precision <- thresholds*0
recall <- thresholds*0

for (k in 1:length(thresholds))
{
  threshold <- thresholds[k]
  cl <- data*0
  tp <- 0
  tn <- 0
  fp <- 0
  fn <- 0

  for (i in 1:length(data))
  {
    if ( (data[i] > threshold) && (i <= Ndata))
    {
      tp <- tp + 1
    }
    if ( (data[i] > threshold) && (i > Ndata))
    {
      fp <- fp + 1
    }
    if ( (data[i] <= threshold) && (i <= Ndata))
    {
      fn <- fn + 1
    }
    if ( (data[i] <= threshold) && (i > Ndata))
    {
      tn <- tn + 1
    }
  }

  precision[k] <- tp / (tp + fp)
  recall[k] <- tp / (tp  + fn)
}

plot(precision,recall)


That really depends on what you need to use the energy for and how exactly you define "energy" in your application.

You can either take the sum of the energies or do a more advanced mode taking into account the inter-channel correlation and the spectral correlation of stereo playback. If you want anything like "physical" energy, you need a calibration.


How to graph magnitude plots for basic 4-pole filters?

For Butterworth LP and HP

https://www.desmos.com/calculator/dxcnvc63av

Play the variable n.

Sources:
https://www.electronics-tutorials.ws/filter/filter_8.html
http://fourier.eng.hmc.edu/e84/lectures/ActiveFilters/node6.html