Difference between Contrast and Intensity of an Image

On intensity

Simply said, it's hard to talk about "intensity" of an image. Every pixel has its intensity (for greyscale images, they are usual allowed range is [0, 255]), but the concept of image intensity does not exist. If you are doing some kind of image analysis, you could be interested in a parameter describing image intensities, e.g. mean intensity (like @geometrikal said) or distribution of image levels (which is related to contrast).

On contrast

If you presume an idealized situation where a greyscale image contains only pixels of two intensity values, one for background pixels and one for object pixels, the contrast of the object would be the difference of those values. What it means on the displayed image --> the higher the contrast, the easier it is to spot (find, locate) the object on the image (for the human eye). Again, you would usually talk about contrast of an object in an image, it's hard to talk about contrast of an image where it is not possible to define an object.

As an example, let's look at several 1-D signals (you can produce images similar to your examples from these):

1)   0   0   0   0  20  20   0   0

2)   0   0   0   0 255 255   0   0

3) 200 200 200 200 255 255 200 200

4) 200 200 200 200   0   0 200 200

I'll assume that the two pixels different from the rest in each example are object pixels, and the rest are background. We have several cases here (and I could make more):

brighter object on dark background, all intensity values are low (i.e. the image is rather dark), contrast is low
bright object on dark background, most pixels have low intensity, contrast is high (it's a black-white image)
bright object on darker background, all intensity values are high (i.e. the image is bright), contrast is low (but higher then in the 1st example)
dark object on bright background, most pixels have high intensity, contrast is high (but not as high as in the 2nd example)

Here's another example supporting the fact that the contrast is usually an attribute describing an object, either in relation to itself or to the surrounding. Another important fact to notice from this example is that contrast generally describes the difference between the intensity levels, but the precise definition depends on the application (purpose for the measurement). Let me paraphrase the contrast definition used in hierarchical image segmentation from P. Soille, L. Najman: On morphological hierarchical representations for image processing and spatial data clustering as an example:

  The internal contrast of a connected component (object) corresponds to the largest intensity difference between two adjacent pixels belonging to this connected component. The external contrast is defined as the smallest intensity difference between a pixel of the considered connected component and an adjacent pixel not belonging to the considered component.

Of course, a different application could use a different measurement for contrast.

On imadjust

I don't usually work in MATLAB, but from the documentation, it is used to map intensity range (intensity levels in given range) of an input image, to an intensity range specified for the output image. The default, no-parameter call, will increase contrast of the image. But, you can get a brighter/darker image with imadjust with either increasing or decreasing the object contrast. Let me demonstrate this on my 3rd example:

imadjust(I3, [0.8; 1.0], [0.0; 1.0])

should output 0 0 0 0 255 255 0 0. You would get a white object on black background. In general, the image is darker (lower intensity), but contrast is higher (black-white).

imadjust(I3, [0.8; 1.0], [0.0; 0.2])

should output 0 0 0 0 51 51 0 0. Fairly dark object on a black background. The image is darker (lower intensity) and the contrast of the object did not change.

imadjust(I3, [0.8; 1.0], [0.85; 0.9])

should output 216 216 216 229 229 216 216. Very bright object on bright background. The image is brighter (higher intensity), and the contrast of the object decreased (object intensity levels are more similar to the intensity levels of the background).

You get the actual ranges the function is working with by multiplying them with the maximum gray level (255). For example, a range [0.2; 0.8] is actually intensity range [51; 204].

One more thing to take care of is that the function clips the values outside the first intensity range, and maps them to the new low if they're smaller or new high if they're larger then the range. All of my examples actually include this: the first range starts from 0.8 which maps to intensity 204, but the intensity of 200 from the input image is mapped to the output low in all the images.

So, it's actually just a simple scaling of image intensities (with cut-offs). Also, the default call to imadjust with only an image as an input parameter should increase contrast. I'd say that imadjust(I2) wouldn't do anything (there's a maximum contrast in my second example).

On contrast enhancement

Quote from P. Soille: Morphological Image Analysis:

  Image contrast enhancement refers to accentuation or sharpening of image features so as to make a graphic display more useful for visualization or analysis of the image by the human eye.

He also emphasizes: with enhanced contrast image analysis by visual (human) inspection is easier. In idealized images like in my examples, a computer wouldn't gain much in term of image analysis difficulty.

E.g. for object extraction by thresholding, it would just mean different thresholds should be used to extract the object from the background. But, a human examiner would spot the object much easier on high-contrast images. This all changes some in real images (objects don't have uniform gray levels), so contrast enhancement becomes useful in image analysis.

There are several methods of contrast enhancement:

point-based techniques, where the local neighborhood is not important. They are based on the analysis of grey levels through the whole image
neighborhood based techniques, where local neighborhood of the pixel is important. Examples are white and black top-hat operators, and toogle contrast operator.
transform based techniques, where the filtering is done on a transformed image before using the inverse transform (e.g. filtering in the frequency domain after Fourier transform).

On intensity adjustment

In term of greyscale image, the intensity of the pixel corresponds to it's brightness. The greater the intensity, the greater the brightness. This also means that increasing intensity can be viewed as brightening the image (while decreasing intensity can be viewed as darkening the image).

I would describe the process of uniformly brightening the image as increasing intensity while leaving the contrast unchanged in the whole image.
This actually means adding a constant value to all the pixels.

Now, as the pixel intensity values have a predefined range, typically [0..255], this means that the maximum intensity exists. This inevitably means that some pixels (that were different from each other before) will become white (i.e. intensity 255).

This "naturally" happens when you take a photo aimed towards something very bright -- sun or other light. Away from the light, you might see some details, but at the position of the light/sun and around it, you will get only white pixels meaning that the intensity (amount of light when the image was taken) "hit" its maximum (displayable/storable) value.

The only way you wouldn't lose details by this operation would be if the original pixel intensities belong to only a part of possible intensity range, e.g. if original pixels are in the range [10, 150], the image could be brightened by up to 105 intensity levels before you start to loose details.

As imgadjust is meant to preform intensity scaling, it can do much more than just brighten/darken image. If you wanted to emulate brightening effect with imadjust, you could write something like:

imadjust(I, [0.0; x], [1.0-x; 1.0])

where x is any number between 0 and 1 (e.g. x=0.5 would brighten a [0..255] image by 127).

That said, this is an overkill for such a simple operation. I'm sure matlab has a elementary operation that adds a scalar to all elements of a matrix, so you could just use that :)


How to Estimate the Input for a Convolution Given the Filter (Impulse Response) and the Output of the Convolution

Mathematically, here is the demostration for what you want to do:
y[n] = x[n]*h[n] (This is a convolution),
Applying convolution theorem you can prove that:
Y[f] = X[f]*H[f] (this is an ordinary product),
Then:
X[f] = Y[f]/H[f],
Now if we apply fourier inverse transform to X[f] you can have x[n].
Implementing this in matlab can be pretty straightforward now.
Cheers


Why are analytic signals so important in Time-frequency analysis?

Assuming time-frequency aims a providing a separation (at least visual) between signal components, the main reasons could be:

for quadratic distributions, which tend to yield interference  between components, "cancelling" at negative frequencies reduce the quantity of components that can interfere.
for linear distributions, the filter bank formalism, and especially the down-sampling operators, is simplified, reducing the impact of aliasing errors.

For real signals, the Hermitian symmetry yield that "no information" is lost in the analytic form, a complex-valued function that has no negative frequency components, and it is easy to go back to real.

However, this is not so simple in practice.

If the analytic signal constructed from a wide-sense stationary (WSS) real signal is always proper, this may not be the case for non-stationary signals, as underlined in Stochastic time-frequency analysis using the analytic signal: why the complementary distribution matters.
Moreover, computing the analytic signal on discrete, finite-length data is often only approximate, especially for real-time applications. Hence the energy in the negative frequencies is rarely zero.
Extending the concept of analyticity beyond 1D is not evident, and several design still exist.


Why do mobile phone calls produce a vibration when they are near a speaker?

A little background...
How a speaker works: Speakers are controlled by driving small diaphrams with a magnet and that magnet is controlled by wires that take isolated 'wired' signals and convert them into sound.

How your cell phone works: Your cell phone has antennas in it.  They work by sending and receiving transmissions of electromagnetic waves. These waves are wireless signals that are picked up by antennas and then processed later.  Your phone is transmitting/receiving those electromagnetic signals.

What is happening between the speaker and the phone: The wires in your speakers are inappropriately acting as antennas and are picking up your phones transmissions.   Note, that it only is happening when the phone is transmitting because it requires significant power to induce a current in your speaker's wires.  Also some consideration must be taken to understand the frequencies that control your speakers but that is more detail than necessary to understand the concept of what is happening.

Read this reference here if you'd like to hear more.


Apply filter before or after Downsampling?

I haven't looked at the reference yet , but just to answer the question, yes filtering is performed prior to down sampling to avoid aliasing back in the out of band signals. However if your signal is already appropriately bandlimited you could forego the filter entirely. There's various reasons you might want to do other, unrelated, filtering after the downsampling operation. In fact, doing such filtering after the downsampler is preferable due to the lower sample rate and thus lower computational complexity.


Normally, in electrical enginnering, we apply the term "transfer function" and "filter" to an operation that belongs in the class we call Linear Time-Invariant systems (LTI).

Sometimes you might read something from someone where they apply either terms to a nonlinear operation.  A 'filter" in statistics usually means a different thing than what EEs mean.  So the use of "filter" is deprecated in this EE context.

Now either a linear or nonlinear operator will take an input signal $x(t)$ or $x[n]$ and will transform it into an output signal $y(t)$ or $y[n]$. The algorithm or rules or math of that operation is what defines what the operator is, linear or nonlinear, time-variant or time-invariant, deterministic or stochastic.

If the operator is LTI, then I would call it either an "LTI system" or a "filter".  For all practical purposes, what defines an operator as linear is

$$\begin{align}
\text{if } \quad y_1(t) &= \operatorname{LTI}\Big\{ x_1(t) \Big\} \\
\text{and } \quad  y_2(t) &= \operatorname{LTI}\Big\{ x_2(t) \Big\} \\
\\
\text{then } \quad  y_1(t)+y_2(t) &= \operatorname{LTI}\Big\{ x_1(t)+x_2(t) \Big\} \\
\end{align}$$

and what defines it as time-invariant is

$$\begin{align}
\text{if } \quad y(t) &= \operatorname{LTI}\Big\{ x(t) \Big\} \\
\\
\text{then } \quad  y(t-\tau) &= \operatorname{LTI}\Big\{ x(t-\tau) \Big\} \quad \text{for any real } \tau \\
\end{align}$$

Now, if your operator satisfies those two conditions, I will call it a "filter".

When a system is LTI, then we can apply all sorts of math on the operator and the signals going in and coming out of the operator.  This includes the Fourier Transform and the Laplace Transform (for discrete-time signals, $x[n]$ and $y[n]$, this would be the Discrete-Time Fourier Transform and the Z Transform, instead).

Probably, the first result we get from Linear System Theory (also called "Signals and Systems" for some textbooks), is that the output $y(t)$ can always be related to the input $x(t)$ by the operation of convolution:

$$\begin{align}
 y(t) &= \int\limits_{-\infty}^{\infty} h(u) x(t-u) \ \mathrm{d}u \\
      &= \int\limits_{-\infty}^{\infty} h(t-u) x(u) \ \mathrm{d}u \\
\end{align}$$

$h(t)$ is the impulse response, and in a sense, fully defines the operation of the filter or LTI system.  I.e. if you know how your LTI system responds to a Dirac impulse as an input, you know how your LTI system will respond to any input.

Then, whether we use the Fourier Transform

$$\begin{align}
 X(f) &= \int\limits_{-\infty}^{\infty} x(t) \, e^{-j 2 \pi f t} \ \mathrm{d}t \end{align}$$

or the (bilateral) Laplace Transform

$$\begin{align}
 X(s) &= \int\limits_{-\infty}^{\infty} x(t) \, e^{-s t} \ \mathrm{d}t \end{align}$$

then the above convolution operation in the time domain becomes a simpler multiplication operation in the frequency domain (or sometimes the "transform domain"):

$$\begin{align}
 Y(f) &= H(f) \, X(f) \\
\\
 Y(s) &= H(s) \, X(s) \\
\end{align}$$

where $X(\cdot)$ and $Y(\cdot)$ and $H(\cdot)$ are all defined the same way from $x(t)$, $y(t)$ and $h(t)$.

The thing that multiplies the (transform of the) input to become the (transform of the) output is the transfer function.  That means $H(f)$ or $H(s)$ transfers the input $X(f)$ or $X(s)$ to the output $Y(f)$ or $Y(s)$.  But that is only in the frequency domain.

In the time domain, the convolution operator (convolving the input $x(t)$ with the impulse response $h(t)$) is what your filter does.

Either the transfer function ($H(f)$ or $H(s)$) or the impulse response ($h(t)$) are sufficient to fully define your LTI system (or filter) from the POV of the input to output relationship.  I.e. two different systems with identical impulse responses $h_1(t) = h_2(t)$ are indistinguishable from each other if both are placed in black boxes.  The same can be said about their transfer functions $H_1(s) = H_2(s)$.  If their transfer functions are identical, from just looking at the inputs and outputs, we cannot tell any difference between the two LTI systems.


In general, an ideal Wiener filter (i.e. a non-causal filter with an infinitely long impulse response) has the following frequency response

$$W(\omega)=\frac{S_{dx}(\omega)}{S_x(\omega)}\tag{1}$$

where $S_{dx}(\omega)$ is the cross-power spectral density of the desired signal $d[n]$ and the input of the Wiener filter $x[n]$, and $S_x(\omega)$ is the power spectral density (PSD) of the input process $x[n]$. If the input is just a noisy version of the desired signal, and if the desired signal and the noise are uncorrelated, then we have

$$S_x(\omega)=S_d(\omega)+S_n(\omega)\tag{2}$$

where $S_n(\omega)$ is the PSD of the additive noise, and $S_d(\omega)$ is the PSD of the desired signal. Furthermore, we have

$$S_{dx}(\omega)=S_d(\omega)\tag{3}$$

because the noise and the desired signal are uncorrelated. Plugging (2) and (3) into (1) gives

$$W(\omega)=\frac{S_d(\omega)}{S_d(\omega)+S_n(\omega)}\tag{4}$$

This expression is very easily interpreted: for frequencies where the noise is much stronger than the signal (i.e. $S_n(\omega)\gg S_d(\omega)$ the Wiener filter strongly attenuates the input signal, whereas for frequencies with a very good signal-to-noise ratio (SNR) (i.e. $S_d(\omega)\gg S_n(\omega)$) the filter passes the input with a gain of approximately $1$.

So your guess is correct: if there's a lot of noise in the high frequencies (relative to the desired signal), then these frequencies are attenuated, and the low frequencies, where the SNR is high, will pass the filter with a gain of approximately $1$.


Finding the gradient of a norm in a minimization problem

This is the result of taking the derivative of $\|\Phi\circ(X_{1}-X_{2})-u\|_F^2$ w.r.t. $X_{1}$:

$[\Phi^{1'}(\Phi\circ(X_{1}-X_{2})-u), ...,\Phi^{J'}(\Phi\circ(X_{1}-X_{2})-u)]$ .

This can be got by first finding the expression for the norm, and then taking the matrix derivative (w.r.t. $X_{1}$) to get an $N \times J$ matrix. This matrix is then represented in terms of the sub matrices $\Phi^{i}$.


Simplified equations for FM (frequency modulation) synthesis and PM (phase modulation)?

If $t$ is time, $s(t)$ is the (appropriately scaled) signal, $\omega_0$ is the angular frequency, and $\phi_0$ is a phase offset, then phase modulation is
$$ t \mapsto \sin(\omega_0 t + \phi_0 + s(t)) $$
Different signs for $s(t)$ and/or $\phi_0$ may be used, depending on conventions and context.

Frequency modulation is more complex to write down, because we want $s(t)$ to vary the derivative of the argument to the sine as a function of time. We get something like:
$$ t \mapsto \sin\bigl(\phi_0 + {\textstyle\int_0^t(\omega_0+s(u))du} \bigr) $$
which is the same as
$$ t \mapsto \sin\bigl(\omega_0t + \phi_0 + {\textstyle\int_0^t s(u)du} \bigr) $$

We see that phase modulation is the same as frequency modulation by the derivative of the signal -- which looks like a boost of high frequencies and attenuation of low frequencies.

Your proposed
$$ t \mapsto \sin\bigl( (\omega_0+s(t))t + \phi_0 \bigr) $$
won't work -- it would look like phase modulation, with a modulation strength that increases monotonically and without bound depending on how long it is since $t=0$.


Not quite the answer you requested, but I think the sanest thing(OOP wise)  to do in most of the situations is to have the backend of the application separated from the frontend.

You could create a standalone library(either static or dynamic) for your dsp code and separate it from the QT code.

This way your DSP code would be independent of the frontend you are using.
You could switch from QT to anything else(wxWidgets/GTK+ etc) and your DSP functions inside the lib would remain the same.


Filtering kinda "half-toned" image for OCR processing

What you really require is probably some morphological operation like dilation followed by erosion. This is called as Closing operation. May be in your case- just dilation itself might be good.

There was a similar question asked previously - which can help with other aspects.

monochrome (1-bit black-and-white) image conversion

How do I reconstruct text from an image using only morphological operations?


Cascading a filter with itself to form a higher-order filter?

Replace your $\alpha$ with the following to see if it helps:

fs = 16000;
fc = 3500;
wc = 2*pi*fc/fs;

L = 5;            % order of the cascade
B = 2^((L-1)/L);
C = cos(wc);
M = B-C-1;
N = 1+C-B*C;
alfa = ( -N + sqrt( N*N - M*M ) )/M;

I've got the following plots, compared to your posted $\alpha$ calculation.


What is the name of this very simple spectral subtraction technique?

Estimating the power of a sum of signals depends on coherence or incoherence assumptions. Details are given for instance on incoherent signal summing or coherent signal summing. Those can provide an estimate of the noise level. They can be called power sums or voltage sums. Other details in Voltage and power add differently. I'll be digging for other names.

The expectation of the energy of a sum $v$ of variables $v_1$ and $v_2$ is:

$$ E((v_1+v_2)^2) = E(v_1^2)+E(v_2^2)+2E(v_1 v_2)\,.$$

Three specific cases. If $v_1=v_2$, you get $+10\log 4 = +6 $ dB. If sources are uncorrelated, $E(v_1 v_2)=0$, so you get  $+10\log 2 = +3 $ dB. If sources are oppositive, they cancel each other, and the total energy becomes $-\infty$ dB.

Then, subtracting the noise level to the observed signal is an instance of (scalar) hard thresholding, called spectral subtraction in the Fourier domain. As eluded to in your Spectrogram with square or non-square magnitude of STFT: power vs. magnitude question, scalar spectral subtraction on the observed STFT ${S}(\tau,\omega)$ can be performed with a power law, as:

$$|\overline{X}(\tau,\omega)|^\alpha=\max\left(|{S}(\tau,\omega)|^\alpha - \lambda|\overline{N}(\tau,\omega)|^\alpha,0\right)\,.$$

The recovery of the noiseless phase is a full different world.


Find signal power in a band by integrating the PSD in a frequency band (and zeroing the out-of-band DFT bins)

For a DFT that was computed with only a rectangular window (no further windowing beyond the selection), the power in the "out of band" bins contains signal energy from "in-band", so by zeroing those out, you are not including the power that should be in the computation. The challenge however is that if there is energy out of band, it will similarly have components within your band of interest! Here we see that the real flaw is trying to comput the power within a bandwidth without windowing the signal in the time domain first to reduce the effects of sidelobes in the DFT.
More details on doing this accurately, and specifically accounting for the power loss due to the windowing, and the change in the equivalent noise bandwidth of each bin, please see this post here:
How can I get the power of a specific frequency band after FFT?
and
How to calculate resolution of DFT with Hamming/Hann window?


Why not use a complex number in the exponent in the z-transform?

$z^{-1}$ is the unit delay operator.

$z$ is the unit advance operator.

$s$ is the derivative operator.

$1/s$ is the integration operator.

the utility of the Z and Laplace transforms is to solve (analyze, predict) constant coefficient differential and difference equations algebraically.


Matrices of complex numbers multiplication

The line
                fscanf(sensorInput, "%lf + j*%lf\t", &x[i]);

should be
                fscanf(sensorInput, "%lf,%lf\t", &x[i].real, &x[i].imag);

Notice that C++ provides complex type, they support the common complex operations, and provide an io, and they can read complex numbers writen as (1,3) directly file_in >> x[i], once you have a istream
Note that this type of question you should ask in stackoverflow instead.


Does 46 dB gain of a filter for a frequency imply 200 times more amplitude?

can we expect to get 200 times higher amplitude in the output of the filter comapre to the input of the filter

Yes.

For a sine wave input the amplitude of the steady state output is simply the amplitude of the input multiplied with the magnitude of the transfer function at that frequency.


A negative frequency is given by the complex signal $e^{-j2 \pi f_c t}$ for which the physical implementation requires two real signal paths, which is made clearer from Euler's formula:
$$e^{-j2 \pi f_c t} = \cos(2 \pi f_c t) + j \sin(2 \pi f_c t)$$
So if the Local Oscillator is generated with sine and cosine outputs, and the multiplier is a complex multiplier (in this case it may have a real input so would result in two real  multipliers) then the receiver can be implemented as described.
Note that a full complex multiplier requires four real multipliers:
$$(I_1+jQ_1)(I_2+jQ_2) = (I_1 I_2 - Q_1 Q_2) +  j(I_1Q_2- I_2Q_1)$$
In this case the input from the LNA is real, representing $I_1$ in the formula above, while the input from the Local Oscillator is complex, representing $I_2+jQ_2$ and the product is:
$$(I_1)(I_2+jQ_2) = (I_1 I_2) +  j(I_1Q_2)$$
With the complex output as indicated in the block diagram.
Note I detail this further and the implications of the real versus imaginary input after the LNA in this post: Frequency shifting of a quadrature mixed signal (for imaginary input the output of the LNA would be fed into a Hilbert Transform and then the multiplier would be a full complex multiplier using four real multipliers).


Bandwidth required to transmit a binary sequence

It's not meaningful, in my opinion, to talk of the "bandwidth" of a constant number.  Unless you want to characterize that constant number as DC (because the value does not change in time).  And the bandwidth of DC is zero.

But a stream of bits toggling on and off, yes, a bandwidth (which is the rate of change of the bits or "bitrate", if we're considering two-sided bandwidth or half the bitrate if we're considering the single-sided bandwidth) is meaningful.

In fact, a very important fundamental theorem from Claude Shannon relates bandwidth, signal-to-noise ratio, and rate of information together.


QAM constellation coder - transmitting a file

Since the bandwidth stays constant at 120 kHz, the symbol rate needs to be reduced: $$R_s = W_m/(1+\alpha) = 120000/2 = 60000 \text{ symbols per second.}$$

Then the bit rate is $R_b = 6R_s=360,000 \text{ bits per second.}$

The transmission duration is then $$T = \frac{48 \times 10^6 \text{ b}}{360 \times 10^3 \text{ b/s}} = 133.3 \text {s}.$$


Get FIR Hilbert coefficients from Matlabs envelope() function

From the mathworks documentation of the function envelope():

  The filter is created by windowing an ideal brick-wall filter with a Kaiser window of length fl and shape parameter $\beta = 8$.

So without hacking the function you can't directly get the filter coefficients, but you can easily find them yourself by just doing what they do, i.e., windowing the ideal impulse response by a Kaiser window. The ideal impulse response of a Hilbert transformer is given by

$$h[n]=\begin{cases}\frac{2}{n\pi},& n\textrm{ odd}\\0,&n\textrm{ even}\end{cases}\tag{1}$$

You just have multiply these coefficients with the window of the desired length.


Obtain a signal's peak value if it's frequency lies between two bin centers

The first algorithm that springs to mind is the Goertzel Algorithm. That algorithm usually assumes that the frequency of interest is an integer multiple of the fundamental frequency. However, this paper applies the (generalized) algorithm to the case you are interested in.

Another problem is that the signal model is incorrect.  It uses 2*%pi*(1:siglen)*(Fc/siglen). It should use 2*%pi*(0:siglen-1)*(Fc/siglen) for the phase to come out correctly.

I also think there is a problem with the frequency Fc=21.3 being very low. Low-frequency real-valued signals tend to exhibit bias when it comes to phase/frequency estimation problems.

I also tried a coarse grid search for the phase estimate, and it gives the same answer as the Goertzel algorithm.

Below is a plot that shows the bias in both estimates (Goertzel:blue, Coarse:red) for two different frequencies: Fc=21.3 (solid) and Fc=210.3 (dashed).  As you can see the bias for the higher frequency is much less.

The plot $x$-axis is the initial phase changing from 0 to $2\pi$.


The upper half of the DFT being the conjugate mirror of the lower half is only true for real valued signals.

For an alternative understanding of what the DFT means you can read my blog article DFT Graphical Interpretation: Centroids of Weighted Roots of Unity.  I also recommend your read my first article "The Exponential Nature of the Complex Unit Circle" for an understanding of the foundation of the DFT, in particular, the significance of the Roots of Unity.

Hope this helps.

Ced


BPSK and QPSK demodulation computational complexity

To demodulate BPSK and QPSK as described by the OP is to simply take the sign of the waveform (decision), which has trivial computational complexity. The real complexity lies in everything that occurs prior to decision in a practical receiver: There is typically an equalizer to remove multi-path distortion, a matched filter to optimize the signal to noise ratio at the decision locations (in the presence of white noise), an automatic gain control (AGC) loop to remove long term amplitude variation, a timing recovery loop to remove phase and frequency offsets in the timing clock (for the proper sample locations of the final decisions), and a carrier recovery loop to remove residual phase and frequency offset of the carrier frequency due to differences in the Tx and Rx clocks and Doppler due to possible motion between transmitter and receiver.
For instance, the block diagram below shows a practical receiver carrier recovery loop that also does the final decisions for QPSK; where for QPSK each decision would be a simple threshold detector on I and Q (notice the relative computational complexity of the decision block versus the rest of the loop, and this is just for carrier recovery and does not include the AGC and timing recovery operations described above that would be typical).


Complex output after inverse FFT of a real signal

Essentially, your code does not respect the inherent Hermitian symmetry of the output of the FFT. Here, your signal is odd-sized $2K+1$. Hence, this FFT yields a complex vector of coefficients $d$ (real) and $a_k$ (generally complex), arranged as:

$$ \left[d,a_1,a_2,\ldots,a_K,\overline{a_K},\ldots,,\overline{a_2},\overline{a_1} \right]$$

If you want to zero-out some coefficients, you should remove then by conjugate pairs (except for $d$), like in

$$ \left[0,0,a_2,\ldots,a_K,\overline{a_K},\ldots,,\overline{a_2},0 \right].$$

Note: for even sizes, a real coefficient $n_{K+1}$ is inserted between $a_K$ and $\overline{a_K}$.

Your line:

X([1:6, 4780:4801]) = 0

does not satisfy this symmetry. Try instead this one (only valid for odd length):

nACBin = 3;
X([1:nACBin+1, end-nACBin+1:end]) = 0; %

and you will get a much lower imaginary residue (here due to numerical issues). Here, it somehow works (see below) because of the very simple context (two sines with lucky interval length), but this is rarely recommended in general:

Especially with short 0 strides, this will convolve your data with a cardinal sine, both broad and wiggling.

As said in comments, you are considering prototyping a filter in the frequency domain using windows. This is a well-covered domain in DSP. However you can avoid those complexities by first looking at filtering in the time domain (Butterworth, Chebyshev, Cauer, Papoulis filters) that already have well-designed implementations. Other filters can be useful for peaks, like polynomial interpolators (Savitzky-Golay, for instance).

Side note: it probably would be more efficient (though a bit more complicated) to embed those frequency constraints into the deconvolution, instead of trying to remove the artifacts afterward.


Determine the haar wavelet function is linear or nonlinear

Linearity can play at different levels in this context:

The Haar wavelet functions are piece-wise affine (constant), not linear (but sometimes people call affine "linear").
Any standard wavelet decomposition is linear, as it decomposes a signal into a sum of coefficients times a wavelet. So this is not specific of the Haar wavelet.
Wavelet approximations are often nonlinear, as they consist in choosing the $K$ best (largest) coefficients to represent the signal. And choosing the $K$ largest does not obey the linearity principle.

So I cannot answer directly, by lack of context.


Hilbert transformer's phase in .AC simulation

You have a linear phase because you probably implemented a causal filter. This just means that you have the desired phase shift of $-\pi/2$ plus a constant delay. The additional delay, which corresponds to a linear phase term, usually causes no problems. It is necessary to make the filter causal and therefore implementable.

The difference between an FIR lowpass filter and a Hilbert transformer is the phase value at frequency zero: it should be zero for a lowpass filter, and $-\pi/2$ for a Hilbert transformer.


Difference between Fourier transform of $\left(\frac{1}{3}\right)^{-k}$ with $k\le -1$ and $\left(\frac{1}{3}\right)^{k}$ with $k\ge 1$

As Atul has pointed out, your both sequences are not the same, but they are time-reversed versions of each other:

k = np.arange(-15, 15)

s1 = (1/3.)**(-k)*(k <= -1);
s2 = (1/3.)**(+k)*(k >= +1)

plt.stem(k, s1, 'r')
plt.stem(k, s2, 'g')

Accordingly, what you see in your calculations is the time-reversal property of the Discrete-Time Fourier Transform:

$$x(t)=y(-t) \Leftrightarrow X(f)=Y(-f)$$

I.e. the frequency changes its sign in your calculations, as is expected from the time-reversal property.


Tradeoff between Filter Order and Sections

Large order single stage filters are inherently problematic. Calculating the numerators and denominator from the poles and zeros is an ill defined problem and numerically unstable.
Furthermore Direct Form II is the worst possible filter topology. I strongly recommend Direct Form I or transposed Form II.

To the questions:

Numerically ill defined problem
Just add the group delay of the indivdiual stages
Group delay is always defined as the negative first derivative of the phase of the transfer function vs. frequency. That's the same for FIR and IIR. It's frequency dependent.
Yes. If you had infinite precision math, the results would be the same.
Very important. You can't implement and unstable filter, you'd only get noise NaN or infinity as output
No idea
Pretty much everyone implements filter ins biquad sections. So the number of stages and the filter order is directly proportional


How to recognize speech bubbles in comic strips (Ideally OpenCV)

As said in the comments an efficient way is to first detect letters, words and text with OCR. Then try to expand each text zone to its corresponding text bubble.

Depending on the text bubble design there are different approaches. However, a solution that could work well and be robust would be to perform edge detection on the near surrounding of the previously detected text zones and then keeping only the (radially) closest edges from the text zones(that should be the bubble borders).

This is just a basic idea, you may have to build a more sophisticated algorithm obviously!


Why there is Difference between shapes of ROC of z domain and s domain?

Because, the region of convergence in the Laplace transform
$$
X(s) = \int_{-\infty}^{\infty} x(t) e^{-st} dt
$$
is related to the weighting provided by the real part of the complex $s = \sigma + j \omega$; as this will yield the weight $|e^{-st}| =  e^{-\sigma t}$  applied on the input signal $x(t)$, and is a function of $\sigma$ alone and is a rectangular (planar) region on the s-plane.

But the region of convergence in the Z-transform
$$
X(z) = \sum_{n=-\infty}^{\infty} x[n] z^{-n}
$$
is related to the weighting provided by the magnitude of the complex $z = \sigma + j \omega$ as given by $|z|^{-n} = |\sigma + j\omega|^{-n} = |z|^{-n} |e^{-j n \angle{z}}| = |z|^{-n}$, which is a circular region on th z-plane, due to the magnitude of $z$ being involved.


Dan Ellis has some neat Matlab scripts that allow you to pull audio features from files . . .

http://labrosa.ee.columbia.edu/matlab/

The dynamic time warp might be a good starting point for the task you describe.


Optimized 2D wavelet transform using FFT

I'm still not convinced you are spending your time optimally. Probably you would get better gains doing other optimizations. Already 10-15 taps are considered very long wavelet filters for most applications that I have seen. But nevermind that I will try and answer your question.

You can do two different approaches

You could do it one level at a time.
Everything at once.

There are drawbacks and gains from either of them. I will start with the first one.

If you go for one level at a time:

You zero pad filters to be equal lenght as the rows and columns. Also make sure you pad so the filter gets in the middle.
FFT the padded filters and corresponding image dimension. If width $\neq$ height you would need to do 2 FFTs per filter one of each dimension length.
Pointwise multiplication of the FFTs of filter and image.
Inverse FFT of the product image.
Downsample product image in the dimension you transformed.
Repeat for the next dimension.
Repeat for all levels.

Notice that this approach will involve a set of iFFTs and FFTs at each level. That would be double-bad.

The everything-at-once approach:

Same as above but instead of iFFT and decimating each level to just FFT back again on the decimated signal you just multiply several times with the filter and do iFFT and decimate 2,4,8 times et-cetera. Then you only need to do one set of FFTs (at the start) but still a full set of iFFTs for each level.

Also the multiplications while in the Fourier domain will be more computationally expensive than ordinary multiplications as the fourier components will be $\in \mathbb{C}$ but filter and image coefficients $\in \mathbb{R}$ ( I guess ) and you would also likely need higher numerical precision to get the same error.


Why multi-path channel has linear phase within the coherence bandwidth?

For the classical baseband channel model for multipath, $h(t)=\sum_ka_k\delta(t-\tau_k)$, the frequency response $H(f)=\sum_k a_ke^{-j2\pi f\tau_k}$. If $K=1$, single tap the phase would have been linear $\phi(f)=2\pi f\tau_k$. But for different path incident on the receiver, there is no guarantee that phase is linear. But let us see how phase changes between first path $k$ and last path $p$. For these incident paths, $2\pi f\tau_k-2\pi f\tau_p \approx \pi$ would cause fading. So if your frequency $f$ changes from $f_0$ to $f_0 +\frac{1}{2(\tau_k-\tau_p)}$, you will encounter fading. Hence $f_w=\frac{1}{2(\tau_k-\tau_p)}$ is called Coherence Bandwidth. It is the amount of bandwidth for which frequency response is 'flat' or constant. If the delay spread $\Delta \tau = (\tau_k-\tau_p)$ was less, $H(f)$ would be varying less over frequencies. If the difference in delay of paths was much much smaller than symbol time, that is, $\Delta \tau << T$, the width of the channel would be much much greater than width of symbols. That is $\frac{1}{\Delta\tau} >> \frac{1}{T}$. So your transmission data would see a 'flat' or constant channel in terms of frequency response. A constant channel in frequency corresponds to a single tap as explained earlier. It's phase response is linear $\phi(f)=2\pi f\tau_k$.


Fourier descriptors: trying to classify objects

The normalization of the Fourier descriptors is performed as follows:

1) Set the DC component of the descriptors to $a(0)=0$.

2) Divide all the Fourier descriptors with the magnitude of the second one. i.e. :

$a(1)=r_1e^{j\phi_1}$

$a(k)=\frac{a(k)}{\lVert a(1) \rVert}$

3) Now, only the position of the starting point remains to be normalized. This is done by subtracting the phase of the second Fourier descriptor $\phi_1$ frin the phase of all Fourier descriptors and weighting by $k$; that is:

$a(k)=a(k)e^{-j \phi_1 k}$

Note that I did not describe rotational invariance here. In simple terms it is achieved by taking the magnitude of each Fourier coefficient.

Also, you have to chain code your contours and represent them in complex space to perform the Fourier transformation (the ordering matters if you don't perform step (2) ).

Fourier descriptors are not necessarily robust. They are good if you have some control of environment (such as lighting and appearance etc. ), but I don't expect them to work on cases such as a missing handle.


Unfortunately, the real-world gets in the way - there's a whole host of things between you and the "bit-perfect" waveform that comes from the MP3 decoder:

As you noted, user controllable filtering (graphic eq for example) - this may be disabled.
"improvements" to the sound made in the digital domain (for example compression to make it sound "louder" and bass boosting) - this may also be disabled
The internal processing is often done at a greater bit depth than the output DAC capability - for example, 32 bit processing but only a 16 bit DAC.  There must then be a process to reduce the bit depth.  This may or may not include dithering of some sort.  The results of this will sound different depending on the implementation.
The digital to analogue converters may be of differing quality, which will impact on (amongst other things) the linearity and noise floor of the analogue output.
The clock feeding the DAC will also have an impact on the performance - higher jitter will be noisier.
The headphone driver will also have various analogue characteristics - eg. noise floor, frequency colouration etc.


Is sub-sampling delay estimation possible without interpolation?

For finite length signals, interpolation (both for frequency and phase) in the frequency domain is helpful if the frequency is not at a DFT/FFT bin center.  Relative phase estimation is required for delay estimation in the frequency domain.  This relative phase estimation can produce a sub-time-domain-sample delay estimation result.  This delay estimation can have an uncertainty related to an unknown integer multiple of the signal's period due to phase unwrapping errors, unless there is some other time domain envelope information available.


Explaining a statement regarding aliasing in a document

The frequencies are $\left|10\,\mathrm{Hz}+k\cdot30\,\mathrm{Hz}\right|, k\in \mathbb{Z}$ and the same for $20\,\mathrm{Hz}$.

Note that real signals always have (complex conjugate) signal components for positive and negative frequencies, so a $10\,\mathrm{Hz}$ signal has components at $\pm10\,\mathrm{Hz}$.  Consequently aliasing in a $30\,\mathrm{Hz}$ grid will also reach $\mp 20\,\mathrm{Hz}$.


I once wrote an application for rectangle detection. It used Sobel edge detection and line Hough transform.

Instead of looking for single peaks in Hough image (lines), the program searched 4 peaks with distance of 90 degrees between them.

For each column in Hough image (corresponding to some angle), three other columns were searched for local maxima. When satifactory peak was found in each of the four columns, the rectangle have been detected.

The program constructed the rectangle and made additional checks for color consistency within and outside the rectangle to discriminate false positives. The program was for detecting paper placement in scanned sheets of papers.


Are reconstruction filter always needed?

If the system naturally has only a low enough frequency response (say a very heavy object relative to the forces involved), and the noise tolerance is high enough, then the effects of any aliasing, due to the lack of any additional low-pass and anti-aliasing filtering, may already be rolled-off below your desired or required noise floor.


What is the Fastest Way Ever to Do 2D Correlation?

In your case, since you have multiple images while you have a given set of kernels the DFT based Correlation would be the best fit.

Pay attention that the DFT Based Convolution / Correlation Is Equivalent to Convolution / Correlation with Periodic / Circular Boundary Conditions.
It means that if you need different boundary conditions (Like padding with Nearest Neighbor) then it means you need to pad the image large enough so the cyclic effect won't hurt you.

In that case you might be better doing, given enough memory, what's being don in DNN.
Apply the Correlation / Convolution with Matrix Multiplication.

Matrix Multiplication is so optimized in our days that even if its complexity is higher, in real world it is faster in many cases.
So you build the Convolution Matrix (You can even chain them) and the images matrix and multiply. I can show a MATLAB Example.

By the way, even if you do it with DFT in the GPU, I suggest you do it with a DNN library (Something which utilize CUDNN, ROCm - Chainer, PyTorch, TensorFlow, MATLAB Neural Network Toolbox, Knet, etc...) since this case of multiple images (Batch) with given kernel is the case optimized in those libraries with automatic mechanism to chose between DFT based or Matrix based implementation.

Good reference - Algorithms for Efficient Computation of Convolution by Karas Pavel and Svoboda David.


Separating waves of very close wavelengths

You have an incorrect approach here. Instead of trying to design small range filters, you should just do frequency-domain analysis.

Say you have a signal $x$ consisting of $N$ samples, sampled at the rate of $f_s$ Hz. For simplicity, let's take an example with $N = 1000$ and $f_s$ = 1000 Hz. The Nyquist frequency is half of the sampling frequency, which is $500$ Hz in this case. This means that the signal you are sampling should not have frequencies above that (otherwise there would be aliasing).

Now, if you want to know what frequencies you have in your signal, compute the discrete Fourier transform (DFT). DFT is defined as

$
X(k) = \sum_{n=0}^{N-1} x(n)e^{-2\pi k\frac{n}{N}} \quad 0 < k < N-1
$

From the DFT, you get two things. First, $|X|$ is called the magnitude spectrum, which is of your interest. The second thing would be the phase spectrum, given by $arg(X)$. DFT is symmetric from the middle for real-valued inputs, so the frequency resolution you get in our example is $1$ Hz. Now, say that you would also want to know the half frequencies. Looking at the definition of DFT, you can see it is $N$-point. If you want, you can do it likewise for $M$-points. Thus, to get the half frequencies as well, take $M=2000$.

In the case of your question, what you need is a really long DFT.

Other things to note: do not compute DFT from its definition. Use a fast Fourier transform (FFT). It's faster, and more precise.


Meaning of the frequencies vector using fft

You forget that what you did was describing the passband signals (i.e. some waveform on a carrier frequency == on a wavelength) in baseband.
So, of course, after calculating the baseband spectrum, you need to shift it back where it came from - to the carrier frequency. It's as easy as that.
If this isn't clear to you, you might want to revisit what (complex equivalent) baseband and passband are.


why is my FFT magnitude of real time recording so volatile?

Are you appending new data to your FFT input vector?  If so, then every FFT result will have a different length; and it's likely the early ones are too short for the frequency resolution you desire.  The smartphone app may also be averaging multiple successive FFT results to reduce noise.


Is there a name for $f(x) = |x|^n \cdot \operatorname{sgn}(x)$?

There are many instances (eg Handbook of Computational and Numerical Methods in Finance, p. 42, Eq. 2.25) where it is simply called the "signed power"  "signed power-$p$ function" or "signed power law".

Some instances are more common. The term "signed square-root" seems to be used, see for instance in this Simulink documentation:

  Square root of the absolute value of the input, multiplied by the sign of the input.

with the signed square as an inverse. It is used in the realm  of statistical distributions, contrast enhancement, or on data features (like SURF) as a compander, close in spirit to the A-law or the mu-law.

For parabolas (and even powers in general), the rewriting into $x\to x.|x|$ could be called a split parabola (see Selected solutions, HW 1, due January 22, 2009, an instance of split functions).


You have to realize that all sequences in DFT relations are inherently periodic. This is usually understood, but it can be made explicit, for instance by that subscript. So both sources are correct, but one explicitly notes the periodicity of the sequence.

If you write down the DFT relation

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

then you see that the expression on the right-hand side of $(1)$ is periodic in $n$ with period $N$. So $\delta[n]$ on the left-hand side is actually continued periodically. If we just consider one period, we don't care about that periodic continuation. However, if we are concerned about the behavior for all $n$, then we do need to take that periodicity into account, and we could have written $\delta_N[n]$ on the left-hand side of $(1)$ to be more explicit.

And in the given example we do care what happens outside that period because the index $l$ in the first sum extends from $-\infty$ to $\infty$.

The last step in the given equation is just a convolution of the sequence $x[n]$ with a unit impulse comb, so $x[n]$ is shifted by integer multiples of $10$ and all those shifted versions are summed up (hoping that the sum converges).


My translation should have described a rotation around the center of the cube, but it was wrong. I changed it to:

rotate45 * (translation - center) + center

and I get the cube in the middle with the right angle.


Conceptual questions from signal processing

Batman has given a great answer. You need to go through the recommended book in order to understand the concepts mentioned. Let me try to simplify it.

BIG PICTURE: De-convolution or inverse filtering is required to retrieve an estimate of the original signal that went through an unknown linear system. Basically, we have a signal which went through an unknown system and we want to get the original back.
Example: You shouted something across a room with lot of reverberation (echos) to your brother. Now for your brother to understand it, he needs to remove the reverberation effect. Let us see it in signal processing terminology.

What is convolution: Let $x(k)$ be your speech signal, and $h(k)$ is the unknown response of the room, then speech heard by your brother is $y(k) = h(k)\ast x(k)$ (this is convolution) and is more formally defined as:

$$y(k) = \sum_{n=0}^{N-1} h(n)x(k-n) $$

Basically, your brother is listening to attenuated copies of your sound, i.e., $y(k) = h(0)x(k) + h(1)x(k-1) + ... + h(N-1)x(k-N+1)$. Note here x(k) is the complete speech of yours and $k-i$ denotes shift by $i$th sample, referring to the shifted copies.

What is deconvolution: Now lets say your brother recorded $y(k)$ and he wants to know what did you say? i.e., what is $x(k)$. What he needs to do is design a filter $w(k)$ such that $\delta(k) = w(k)\ast h(k)$, i.e., $w(k)$ is the inverse of $h(k)$. If you know basic signal processing then $\delta$ is the dirac delta function or an impulse. Now if we convolve $y(k)$ with $w(k)$ then $w(k)\ast y(k) = w(k)\ast h(k) \ast x(k) = \delta(k) \ast x(k) = x(k)$. And he can know what you said.

Minimum Phase: deals with the design of $w(k)$. Minimum phase also means that our filter is designed such that it responds much faster than another filter that is non-minimum phase. It has to do with setting of zeros of $w(k)$. This concept will require another discussion, which I'll postpone for now.


The impulse invariance method as well as the bilinear transform will give you relatively large errors close the Nyquist and there is nothing you can do about it. One thing you could do is use a filter design method that approximates the given analog response. I've tried this for the RIAA de-emphasis fitler with an (unpublished) algorithm I came up with some time ago, and the figure below shows the result:

The digital filter is a second-order IIR filter, the sampling frequency is $f_s=44100\text{ Hz}$, and the error is much less than 1dB over the whole frequency range up to Nyquist. The price you pay for it is 20 iterations of an approximation method, which solves a $5\times 5$ system of linear equations in each iteration.

EDIT:

If you add one more pole and one more zero (i.e. your discrete-time filter has 3 poles and 3 zeros), then the error goes down to less than 0.1 dB over the whole frequency range up to Nyquist:


You can find the formula for (at least) one "Harris method" or "Harris approximation" in Multirate Signal Processing for Communication Systems, Fredric J. Harris, 2004, page 216, equation  (8.16), with details in Chapter 3. F. Harris calls it "harris approximation" (without capital for the name), namely:

$$ N = \frac{f_{SMPL}}{\Delta f}\frac{Atten(dB)}{22}$$

Also in Digital Communications with Emphasis on Data Modems: Theory, Analysis, Richard W. Middlestead, page 752, equation B.39.

Former questions might provide you with additional insights and references:

How many taps does an FIR filter need?
Filter Order Rule of Thumb


Can one guess what stage of experiment every matrix refers to? If so, how? I mean, do they at least look like the result of the experiment described? If they do, how can one interpret them?

I do not think that I can. But here is a guess that might give you some insight: The most "populated" matrices are those for a single hand up in the air (whether left or right).

Why?

Functional connectivity metrics (of which there are so many by now) assess some of "dependence" or "correlation". There is a higher chance for any two pairs of waveforms to have more correlated signals when the diversity of the generated signals is low. With two hands being raised, the brain generates a large set of impulses to drive all the necessary muscles. Consequently, the chance of them having a high correlation is low.

If you notice, a key characteristic of functional connectivity estimation in related papers is some form of averaging, whether temporal, spatial or across and within subjects.

Whether it is Magnetoencephalography (higher spatial resolution than Electroencephalography), or of course, fMRI or EEG, there is an intermediate stage that "summarises" observations.

In fMRI you sum voxel traces around anatomical areas indicated by an atlas, in MEG you sum contributions of elements around broad areas (simplest method) OR you do some form of source localisation and then look at the connectivity of the localised sources. Similar workflow applies in EEG too (e.g. LORETTA as above).

This is a required step to increase the Signal to Noise Ratio (SNR) and look at (hopefully) robust correlations. Along this idea are also methods that examine functional connectivity in specific bands of brain activity

Therefore, I doubt that you will get anything meaningful by looking at each individual channel without some form of spatial or temporal filtering / averaging.

This however, is half of the story. The other thing you can do is, of course, to change the task.

  What could be an experiment that would show the noticeable PSI for two areas given beforehand? I mean, we had two areas given, pre-motor cortex and primary motor cortex, and we supposed that we would find a synchronization between them in such a motion experiment; something's gone wrong, though. If the problem is in the design of experiment, what design would be better?

You will definitely see differences in the default mode network between two radically different tasks. For example, at rest and at "work".

To reduce your post-processing workload, try to use tasks that are not likely to generate a lot of artifacts. For example: Resting state and "solve a mathematical problem".

This design will give you something to compare "activity" against and help in reducing spurious correlations further.

Finally, be aware of filtering. The more your filter a waveform, the more it "looks like the filter" and therefore, the more chances there are to have more increased correlations.

Hope this helps.


Is ideal band pass filter (brick wall filter) linear phase?

Assuming you are refering to LTI (linear time-invariant) systems to implement the filter.

The impulse response $h[n]$ of the ideal brickwall bandpass filter :
$$H(\omega) = \begin{cases} 1 ~~~,~~~ |\omega-\omega_c|<W \\ 0 ~~~,~~~\text{o.w.}\\ \end{cases}$$
is
$$
h[n] = 2 \cos(\omega_c n) \frac{ \sin(W n) }{ \pi n }
$$

which is real and even symetric about origin; i.e., $h[n] = h[-n]$,  hence it's a zero-phase (in the passband), but noncausal filter.

When this impulse response is truncated and shifted right enough to make it causal: $h_c[n] = h[n-d]$ then the resulting bandpass filter will be linear phase in the passband:
$$ H_c(\omega)= e^{-j\omega d} H(\omega)$$ with a linear phase term of $$\phi(\omega) = -\omega d$$

The phase is not defined for those frequencies where the frequency response is zero.

Note that the truncated filter is no more the ideal filter. Also note that you cannot shift the ideal filter to make it causal (requires infinite shift). And finally note that the zero-phase ideal filter is also a linear phase filter.


Answers to Your Questions

You are right that your time instants are spaced at an interval of $T$. Hence, the sampling rate is $1/T$. ($T$ is defined in line 3 of your code.) You are also right that the second parameter of np.fft.fftfreq() is the sampling interval $T$. I am, however, not sure why you get wrong results when setting $T=100\,\text{Hz}$.
This phenomenon is called aliasing. Nyquist's sampling theorem states that frequencies below $f_\text{s}/2=\frac{1}{2T}$ can be captured correctly by the sampling process. Frequencies above half of the sampling frequency are mirrored back to frequencies below half of the frequency. Therefore, a frequency of $1080\text{Hz}$ is too high to be captured correctly.
The choice for the second parameter of the FFT, the FFT length, depends on the specific application. You have to experiment what leads to the best results. Still, your choice of NFFT = 1024 is typically a good start.


How to discover the caracteristic frequency of the use of an app by someone

I think a Fourier based frequency spectrum analysis may not be the best (most suitable) approach in providing you the necessary and useful information about that random service usage data.

Instead I would suggest you to use statistical data analysis methods. These methods are the primary tools that are used to analyse user behaviour by collecting usage data and computing the resulting statistics upon them.

Creating the related random variables and estimating their associated probability density functions will allow you to formulate the average quantities such as the expected number of days (or hours) between two consequetive logins, expected day of week (or hour of day) of maximum usage, minimum usage etc, and therefore enable you to make optimizations of resource utilizations based upon them.


Why delete high frequence in fft return the area of the most changes?

This line of yours

  delete the component with magnatude compare to DC component(=2.5e5),
  i.e. if Af(i,j)>1e3 , then Af(i,j)=0

indicates an effective highpass filter on the typical smooth image data.

Indeed for typical smooth images high, frequency coefficients are generally (but not always) smaller in magnitude compared to the low frequeny coefficients. So deleting large magnitude coefficients will effectively delete low frequency coefficients. But this is not exactly same thing with deleting low frequency coefficients directly. neverthless, such an operation would return mostly a darkened image with certain edges and contours outlined, essentially a very a primitive edge detection exemplified with the following code:

I = im2double(imread('Lena.bmp'));
figure,imshow(I);

Ik = fft2(I);
th = 100;
Ik( abs(Ik) > th) = 0;
figure,imshow( 3*real(ifft2(Ik)));

with the result:


The frequency response of the two complementary filters are $H_2(e^{j\theta}) = 1 - H_1(e^{j\theta})$, or the impulse responses $h_2[n] = \delta[n] - h_1[n]$.

For an IIR filter, $H_1(z)$ can be written as $\frac{b_0 + b_1 z^{-1} + \ldots}{a_0 + a_1 z^{-1} + \ldots}$. Then $H_2(z)$ should be something like $\frac{(a_0 - b_0) + (a_1 - b_1) z^{-1} + \ldots}{a_0 + a_1 z^{-1} + \ldots}$.

So the non-recursive coefficients for $H_2$ are now $(a_0-b_0)$, $(a_1 - b_1)$, etc.

The recursive coefficients are the same for both filters.


autoregressive moving average code implementation

I dont know in Python, but this simple matlab/octave loop code:

function y=arma(x,A,B)
na=length(A);
n=length(x);
y(1:na,1)=zeros(na,1);
for i=na:n
    y(i,1)=1/A(1)*A(2:na)*y(i-1:-1:i-na+1)+1/A(1)*B*x(i:-1:i-na+1);
end

Implements a ARMA filter, with polynomials $A$, $B$, such as per this 5th order filter:

A=[1 0 0 0 0 0];
B=[0 0.2 0.2 0.2 0.2 0.2];
n=1000; dt=0.001;
t=dt*(0:n-1)';
x=randn(n,1);
y=arma(x,A,B);
plot(t,[x y]);


Your description of how to create a spectrogram is correct. Here's an explanation with a little more detail:

First, create a WaveChannel32, using a WaveFileReader. This will output a PCM byte stream. Then read 4 times the size of your FFT bytes from the WaveChannel32 into a buffer. Convert every four bits in the buffer into an integer.

It's important at this point to use a window function on our interval, to emphasise amplitudes closer to the current time stamp.

I'd then suggest using the Exocortex FFT library to calculate the FFT. The Fourier.FFT method will convert an array of complex numbers to its corresponding Fourier transform. Then ignore the second half of the array produced, to prevent aliasing.

The first half of the array corresponds to the column of the spectrogram at our current time stamp. Each number in the array corresponds to a single pixel in the spectrogram, the frequency of which is determined by the index of the array (equally spaced throughout 0 to 22050Hz). The darkness of each pixel is determined by the magnitude squared of the complex number corresponding to said pixel.

Here's some sample code for producing one column:

var waveChannel = new NAudio.Wave.WaveChannel32(new NAudio.Wave.WaveFileReader("file path"));

var buffer = new byte[4096*4];

waveChannel.Read(buffer, 0, 4096*4);

var samples = new ComplexF[4096];

for (var i = 0; i < 4096, i++)
{
    samples[i].Re = BitConverter.ToSingle(buffer,i*4);

    // Window should be pre-populated with values from your window function
    samples[i].Re *= Window[i];
}

Fourier.FFT(samples, FourierDirection.Forward);


IIR coefficients and difference equation implementation in C language

/**********************************************************************
*                                                                     *
*                                                                     *
*                                                                     *
*                        Copyright (c) 2005                           *
*                                                                     *
*                      Robert Bristow-Johnson                         *
*                                                                     *
*                     rbj@audioimagination.com                        *
*                                                                     *
*                                                                     *
*                                 EQ.c                                *
*                                                                     *
*                                                                     *
*   Utilities for computing EQ filter coefficients, plotting          *
*   frequency response, and EQ filtering sampled audio data.          *
*                                                                     *
*                                                                     *
**********************************************************************/

#include <stdio.h>
#include "__functions.h"

enum EQ_type
{
LPF,
HPF,
BPF,
BPF0,
notch,
APF,
peakingEQ,
lowShelf,
highShelf,
lowShelf1,
highShelf1
};

typedef struct
{
float Fs;
enum EQ_type type;
float f0;
float dBgain;
float Q;
float filterCoefficients[5];    /* b0, b1, b2, a1, a2   (a0 is normallized to 1) */
float graphCoefficients[6];     /* B2, B1, B0, A2, A1, A0 */
} sectionData;

int main(void);
int computeSectionCofficients(float Fs, enum EQ_type type, float f0, float dBgain, float Q, sectionData *data);
int clearArray(long N, float *array);
int plotSectionFrequencyResponse(int Npoints, float Fmin, float Fmax, float *frequencyValues, float *dBvalues, sectionData *data);
int filterBuffer(long Nsamples, int Nsections, float *input, float *output, float *filterStates, sectionData **sections);

#define BUFFER_SIZE 1024
#define GRAPH_SIZE 589
#define NUM_SECTIONS 4
#define FS 44100.0
#define FMIN 19.99955848414893
#define FMAX 19999.55848414893

#define F1 50.0
#define G1 0.0
#define Q1 1.0

#define F2 200.0
#define G2 0.0
#define Q2 2.0

#define F3 1000.0
#define G3 10.0
#define Q3 2.0

#define F4 5000.0
#define G4 0.0
#define Q4 1.0

sectionData section1, section2, section3, section4;         /* this has to be global */

void plotTheFrequencyResponse(void);
void runTheFilter(void);

int main(void)
{
computeSectionCofficients(FS, lowShelf1,  F1, G1, Q1, &section1);
computeSectionCofficients(FS, peakingEQ,  F2, G2, Q2, &section2);
computeSectionCofficients(FS, peakingEQ,  F3, G3, Q3, &section3);
computeSectionCofficients(FS, highShelf1, F4, G4, Q4, &section4);

plotTheFrequencyResponse();
runTheFilter();
}

void plotTheFrequencyResponse(void)
{
float freq[GRAPH_SIZE];
float graph[GRAPH_SIZE];

clearArray(GRAPH_SIZE, graph);
plotSectionFrequencyResponse(GRAPH_SIZE, FMIN, FMAX, freq, graph, &section1);
plotSectionFrequencyResponse(GRAPH_SIZE, FMIN, FMAX, freq, graph, &section2);
plotSectionFrequencyResponse(GRAPH_SIZE, FMIN, FMAX, freq, graph, &section3);
plotSectionFrequencyResponse(GRAPH_SIZE, FMIN, FMAX, freq, graph, &section4);
}

void runTheFilter(void)
{
float aFilterBuffer[BUFFER_SIZE];
float aFilterStates[(2*NUM_SECTIONS+2)];
sectionData *sectionPtrArray[NUM_SECTIONS];

sectionPtrArray[0] = &section1;
sectionPtrArray[1] = &section2;
sectionPtrArray[2] = &section3;
sectionPtrArray[3] = &section4;

clearArray(2*NUM_SECTIONS+2, aFilterStates);
clearArray(BUFFER_SIZE, aFilterBuffer);
aFilterBuffer[0] = 1.0;                         /* input an impulse to get impulse response */
filterBuffer(BUFFER_SIZE, NUM_SECTIONS, aFilterBuffer, aFilterBuffer, aFilterStates, sectionPtrArray);
}

int computeSectionCofficients(float Fs, enum EQ_type type, float f0, float dBgain, float Q, sectionData *data)
{
register float temp, *coef_ptr;
float A = __exp2(0.08304820237218*dBgain);
float A2 = A*A;
float Ap1, Am1;
float cos_w0 = __cos_pi(2.0*f0/Fs);
float sin_w0 = __sin_pi(2.0*f0/Fs);
float alpha = 0.5*sin_w0/Q;
float a0, a1, a2, b0, b1, b2;

switch (type)
    {
    case LPF:
        b0 =  0.5*(1.0 - cos_w0)*A2;
        b1 =  2.0*b0;
        b2 =      b0;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case HPF:
        b0 =  0.5*(1.0 + cos_w0)*A2;
        b1 = -2.0*b0;
        b2 =      b0;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case BPF:
        b0 =   0.5*sin_w0*A2;
        b1 =   0.0;
        b2 =  -b0;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case BPF0:
        b0 =   alpha*A2;
        b1 =   0.0;
        b2 =  -b0;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case notch:
        b0 =   A2;
        b1 =  -2.0*cos_w0*A2;
        b2 =   A2;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case APF:
        b0 =  (1.0 - alpha)*A2;
        b1 =  -2.0*cos_w0*A2;
        b2 =  (1.0 + alpha)*A2;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case peakingEQ:
        if (A > 1.0)
            {
            alpha *= A;
            }
         else
            {
            alpha /= A;
            }
        b0 =   1.0 + alpha*A;
        b1 =  -2.0*cos_w0;
        b2 =   1.0 - alpha*A;
        a0 =   1.0 + alpha/A;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha/A;
        break;

    case lowShelf:
        temp = 2.0*__sqrt(A)*alpha;
        Ap1  = A + 1.0;
        Am1  = A - 1.0;
        b0 =  A*(Ap1 - Am1*cos_w0 + temp);
        b1 =  2.0*A*(Am1 - Ap1*cos_w0);
        b2 =  A*(Ap1 - Am1*cos_w0 - temp);
        a0 =  Ap1 + Am1*cos_w0 + temp;
        a1 =  -2.0*( Am1 + Ap1*cos_w0);
        a2 =  Ap1 + Am1*cos_w0 - temp;
        break;

    case highShelf:
        temp = 2.0*__sqrt(A)*alpha;
        Ap1  = A + 1.0;
        Am1  = A - 1.0;
        b0 =  A*(Ap1 + Am1*cos_w0 + temp);
        b1 = -2.0*A*(Am1 + Ap1*cos_w0);
        b2 =  A*(Ap1 + Am1*cos_w0 - temp);
        a0 =  Ap1 - Am1*cos_w0 + temp;
        a1 =  2.0*( Am1 - Ap1*cos_w0);
        a2 =  Ap1 - Am1*cos_w0 - temp;
        break;

    case lowShelf1:
        temp = cos_w0 + 1.0;
        b0 =  A2*sin_w0 +    temp;
        b1 =  A2*sin_w0 -    temp;
        b2 =  0.0;
        a0 =     sin_w0 +    temp;
        a1 =     sin_w0 -    temp;
        a2 =  0.0;
        break;

    case highShelf1:
        temp = cos_w0 + 1.0;
        b0 =     sin_w0 + A2*temp;
        b1 =     sin_w0 - A2*temp;
        b2 =  0.0;
        a0 =     sin_w0 +    temp;
        a1 =     sin_w0 -    temp;
        a2 =  0.0;
        break;
    }

coef_ptr = &(data->filterCoefficients[0]);
temp = 1/a0;
*coef_ptr++ = temp*b0;                      /* b0/a0 */
*coef_ptr++ = temp*b1;                      /* b1/a0 */
*coef_ptr++ = temp*b2;                      /* b2/a0 */
*coef_ptr++ = temp*a1;                      /* a1/a0 */
*coef_ptr   = temp*a2;                      /* a2/a0 */

coef_ptr = &(data->graphCoefficients[0]);
temp = b0*b2;
*coef_ptr++ = temp;                         /* B2 */
*coef_ptr++ = -0.25*b1*(b0 + b2) - temp;    /* B1 */
temp = b0 + b1 + b2;
temp = 0.0625*temp*temp;
*coef_ptr++ = temp;                         /* B0 */
temp = a0*a2;
*coef_ptr++ = temp;                         /* A2 */
*coef_ptr++ = -0.25*a1*(a0 + a2) - temp;    /* A1 */
temp = a0 + a1 + a2;
temp = 0.0625*temp*temp;
*coef_ptr   = temp;                         /* A0 */

data->Fs = Fs;
data->type = type;
data->f0 = f0;
data->dBgain = dBgain;
data->Q = Q;

return 0;
}

int clearArray(long N, float *array)
{
register long n;

for (n=0; n<N; n++)
    {
    array[n] = 0.0;
    }

return 0;
}

int plotSectionFrequencyResponse(int Npoints, float Fmin, float Fmax, float *frequencyValues, float *dBvalues, sectionData *data)
{
register int n;
float logFreqStep = __log2(Fmax/Fmin)/(float)(Npoints-1);
register float phi, tau = 1.0/(data->Fs);
float *coef_ptr = &(data->graphCoefficients[0]);
register float B2, B1, B0, A2, A1, A0;

B2 = *coef_ptr++;
B1 = *coef_ptr++;
B0 = *coef_ptr++;
A2 = *coef_ptr++;
A1 = *coef_ptr++;
A0 = *coef_ptr;

for (n=0; n<Npoints; n++)
    {
    phi = Fmin*__exp2((float)n*logFreqStep);
    frequencyValues[n] = phi;
    phi = __sin_pi(tau*phi);
    phi = phi*phi;                                                                              /* sin(w/2)^2 */
    dBvalues[n] += DB_LOG2_ENERGY*(__log2((B2*phi + B1)*phi + B0) - __log2((A2*phi + A1)*phi + A0));
    }

return 0;
}

int filterBuffer(long Nsamples, int Nsections, float *input, float *output, float *filterStates, sectionData **sections)
{
register long n;
register int i;
register sectionData **section_ptr;
register float sampleValue, state1, state2, *state_ptr, *coef_ptr;

for (n=0; n<Nsamples; n++)
    {
    section_ptr = sections;                         /* reset section pointer */
    state_ptr = filterStates;                       /* reset state pointer */

    sampleValue = input[n];                         /* get input sample */

        state2 = *state_ptr++;                      /*   x[n-2]            */
        state1 = *state_ptr--;                      /*   x[n-1]            */
        *state_ptr++ = state1;                      /*   x[n-1] -> x[n-2]  */
        *state_ptr++ = sampleValue;                 /*   x[n]   -> x[n-1]  */

    i = Nsections;
    while (--i >= 0)
        {
        coef_ptr = &((*section_ptr++)->filterCoefficients[0]);  /* point to section filter coefficients */

        sampleValue =  *coef_ptr++ * sampleValue;   /*   b0*x[n]           */
        sampleValue += *coef_ptr++ * state1;        /*   b1*x[n-1]         */
        sampleValue += *coef_ptr++ * state2;        /*   b2*x[n-2]         */

        state2 = *state_ptr++;                      /*   y[n-2]            */
        state1 = *state_ptr--;                      /*   y[n-1]            */
        *state_ptr++ = state1;                      /*   y[n-1] -> y[n-2]  */

        sampleValue -= *coef_ptr++ * state1;        /*  -a1*y[n-1]         */
        sampleValue -= *coef_ptr++ * state2;        /*  -a2*y[n-2]         */

        *state_ptr++ = sampleValue;                 /*   y[n]   -> y[n-1]  */
/*
*      at this point sampleValue is the output of the section just completed
*     and will be the input to the section about to be computed unless this
*     is the last section, in which case sampleValue is the input.
*     state1 and state2 are the output states of the section just completed
*     and will be the input states of of the section about to be computed.
*/
        }

    output[n] = sampleValue;                        /* store output sample  */
    }

return 0;
}


Sampling from a distribution with arithmetic coding

While this answer may have come a bit too late for OP, here it is.

Let's start by replacing the phrase 'bits used' with 'queries to an unbiased bit generator'. Then assume that the only random number generator you have at your disposal is an unbiased bit generator, that generates $0$ or $1$ with probabilities $p_0=p_1=0.5$. e.g. something like the following (well, this is probably not the best way to implement an unbiased bit generator in octave/matlab, but suffices for our purpose):

function b = unbiased_bit_generator() b=randi([0, 1], 1,1); end;

Now your goal is to use this function as less as possible (1 use = 1 query) to generate a given number of iid biased bits with a specified bias.

In case of rejection sampling, you would use $\texttt{unbiased_bit_generator()}$ 32 times to get a 32-bit integer $i$, then check if $\frac{i}{2^{32}}<p_0$ to generate one biased bit. Thus you require $32$ queries/biased-bit. Here is a biased bit generator based on 32-bit rejection sampling:

function b=biased_bit_generator1(p0)
    for i=1:32
        bin_array(i)=unbiased_bit_generator();
    end;
    str=num2str(bin_array);
    str=str(str~=' ');
    if bin2dec(str)/2^32<p0
        b=0;
    else
        b=1;
    end;
end;

In case of arithmetic decoder based sampling, you feed the decoder with a sequence $s$ of $k$ unbiased bits (generated using $\texttt{unbiased_bit_generator()}$ $k$ times) and get a sequence $d$ of of length $n$ biased bits. Thus you needed $\frac{k}{n}$ queries/biased-bit.  Because arithmetic coding is a variable-length coding technique, two encoded sequences of same length may correspond to two message sequences of different length. Thus, if you were to end up with  sequence $s_k'\neq s_k$ after using the unbiased bit generator $k$ times and feed it to the arithmetic decoder, then you may get as output a sequence $c_{n'}$ s.t. $n\neq n'$.  In this case, you needed $\frac{k}{n'}$ queries/biased-bit.

The number 0.081 specified in the book is the expected number of queries/biased-bit where the expectation is over all encoded input sequences to the decoder. As you keep feeding your decoder with more and more unbiased bits and keep generating more and more biased bits the queries/biased-bit will approach this number.

  Also this seems a bit too good to be true

I assume you are talking about 0.081 queries/biased-bit. I think this statement is equivalent to saying that Shannon's source coding theorem is too good to be true ;)


I think you misunderstand what this filter is doing. It does not compute approximations of the Bessel functions $J_0(x)$ and $J_1(x)$, but it computes an approximation of the Hankel transform of some input function. The function is expressed as a weighted sum of Bessel functions, and the result of the transform are the weights. The digital filtering approach is an efficient alternative to solving the transform integral by standard methods of numerical integration.


Would it be possible to reconstruct 3d trajectory of ball projectile using only monocular video

I guess this is what Olli implies: Using a lens distortion matrix if necessary, the center of the ball in the picture gives you an elevation and azimuth of your ball relative to the focal point of your camera. The size of the ball image gives you distance.

That's a point in 3D space.

From there, it's really only a manner of grammar school level math to estimate the four (or only three – do you know the mass of the planet you're on?) unknowns in your problem (three components of initial velocity, and maybe magnitude of gravity), at least without considerations for measurement error from multiple observations of the ball.

If you need to take measurement error (finite shutter speed, finite frame rate, finite camera resolution, imperfect estimation of position from picture) into account, well, feels like the classical problem where you'd throw a Kalman at it and be done.

In fact, I think this nice little article "How a Kalman filter works, in pictures" has an explanation that nearly completely describes a solution for your specific problem.


Relative merits of fixed point vs floating point computation?

The numerical accuracy of integers will only be better than the numerical accuracy of floats if the integer resolution is better.  Doubles have 52 fractional bits, so double-precision floats have a resolution worse than integers at around $2^{52}$, which is much larger than 32768 ($2^{15}$).  So, no, the numerical accuracy will not be better if you go to integers.

The second issue is speed.  The answer is: it depends on the hardware.  If you are running your program on a digital signal processor that has multiple fixed point multiply/accumulate cores then yes, it will be far faster in fixed point.  On an x86 chip, on the other hand, it will actually probably be slower in fixed point.  I did exactly what you are talking about once and saw my execution times increase.

After doing some searches on the internet I found out that that is common.  The reason is because it has a dedicated floating-point processor that isn't doing anything when you transition to fixed point, while the fixed-point hardware is shared with the regular fixed-point action, such as pointer arithmetic.

If you want to speed up the processing, the way to do it is to change from double-precision floats to single-precision floats.  That should produce a significant increase in speed.  That would, of course, reduce your numerical accuracy though.


What you'd often be looking for would the variance – the expectation of the squared difference to the mean value in that region. If the noise you add to the actual image is zero-mean, then the mean of the flat region is (in expectation) the actual intensity, and subtracting it from all pixels and squaring the result would give you the noise power.

However, shot noise is not zero-mean:

each flat-zone image pixel $v_i$  follows an actual-intensity-"offset" poisson distribution

$$v_i = m + n_i,\quad n_i\sim\text{Pois}(\lambda)$$

with $m$ being the actual intensity of the object, $n_i$ being shot noise samples, and $\lambda$ the intensity of the Poisson variable.

Now, since shot noise and actual image are independent,

$$\mathbb E(v_i) = \mathbb E(m + n_i) = \mathbb E(m) +  \mathbb E(n_i) =   m + \lambda\text.$$

Sadly, since we don't know the actual $m$  a priori, we'll have to dig one moment  deeper:

\begin{align}
\DeclareMathOperator{\Var}{Var} \Var (v_i) &= \mathbb E\left((m + n_i- \mathbb E(v_i))^2\right)\\
& = \mathbb E\left((m + n_i- m - \lambda)^2\right) \\
& = \mathbb E\left(( n_i - \lambda)^2\right) \\
&=: \Var(n_i)\\
&=\lambda
\end{align}

So, with the Variance of your observation sample as estimator for the $\lambda$ of your Poisson variable, you get the sole statistical property of your shot noise for free. If you want to know the actual flat color, you'd also subtract that from the average of the flat region.


Does this get you anywhere? Apologies for my lacking coding vocabulary, feel free to edit.

Figure 1. Adobe Audition spectrogram of the lowest band of 06-Sep-2018 18h48m47s - 460.519.5 MHz, NFM.wav using Kaiser (180dB) window, 2048 bands, 100 % window width, 80 dB range. A manually adjusted grid is superimposed over the spectrogram.

It's repeating a 16-bit pattern ???? ???? ?111 1110 (spaces added for readability). That would be 9 data bits, and a synchronization sequence of 6 up bits and 1 down bit. The data bits should not contain the string 1111110 or else synchronization may fail. You can synchronize the bit clock to the first up→down transition after a long enough run of the up frequency.


How to detect peaks ignoring soft declivity?

You mentioned that you have tried the signal derivative to catch the declines. What problems have you faced?

It seems that your problem might be addressed with different words: Let $s(t)$ be your signal, and let $u(t)=ds(t)/dt$. You want a detector to indicate when $a\le u(t)\le b$, with $a,b<0$, for at least $T$ seconds.

If this rewrite of the problem is true, the challenge now is to calculate $u(t)$. You can use u=diff(s), as you mentioned if the signal $s(t)$ has no noise, or you need to estimate $u(t)$ using other techniques like these two for instance.


Adaptive equalizer for BPSK signals with frequency offset

It sounds like your biggest issue is dealing with a phase offset and otherwise you have no other concerns with your equalization approach?  In that case, the issue is resolved by simply having a phase correction prior to making your decisions on the real signal. If done correctly your carrier recovery implementation should remove any residual long term phase offset. This post may help give you further insight: Phase synchronization in BPSK.  I duplicated the final block diagram below where in your case "I data demod" would be your BPSK output. I want to point out using this diagram, that any residual phase error would accumulate in the PI Loop Filter which would then adjust the NCO to rotate the signal to have zero error (on average). Note that the Loop as shown if completely blind will lock in quadrature or 180° phase positions - if that is your remaining issue then another question on carrier recovery implementations would be more appropriate specific to any issue you have there, but a simple solution is to derive the phase error itself from $\phi \propto Q \text{sign}(I)$ since the objective in this case is to drive any signal on Q down to the noise. The assumption is the signal has been AGC'd prior to this point so such an approach to phase detection would result in a consistent loop gain.  Here is another related post with regards to what gain constants you should use (what loop BW should be used): Loop bandwidth for symbol timing recovery

Also notice the similarity of the above diagram to a classical Costas Loop with the following figure from http://michaelgellis.tripod.com/mixerscom.html. This is because the Decision Directed CR Loop is a Costas Loop. The referenced link also diagrams the classical Costas Loop for BPSK specifically, which also shows how the digital implementation such as shown in the diagram above can be modified for the BPSK only case (which is modified the phase error metric as I described earlier).  This solution will still have a 0/180° phase ambiguity which is typically resolved with minimum knowledge of the message (header or some other recognized feature in the data message or encoding in the data message).


How is binary data related to double precision data?

I don't know much about Matlab's wavread, but double precision is a standard for describing a "floating point" value. That is, a value that always has the same number of (binary) digits of precision.
A double precision floating point number in binary looks like this:
Each of those boxes represents a bit (short for "binary digit", or a 1 or a 0). It works like scientific notation: the exponent and fraction are both base-2 representations of numbers, and the number this represents is (with some exceptions for things like 0 and NaN) $$(-1)^{\text{sign}}\times 2^{\text{exponent - exponent bias}} \times 1.\text{fraction}$$
Where the "exponent bias" is there to allow you to represent negative exponents, and the "." in $1.\text{fraction}$ is a "binary point" instead of a decimal point.
You can read about the particulars of the double precision format here, which is where I got the image and formula.


Linear model of the modified error function

Yes, you can do that.

$$E(w)=1/2\sum_{n=1}^{N}\{y(x_n,w)-t_n\}^2+\dfrac \lambda2||w||^2$$

$$\begin{align} \frac{\partial E(w)}{\partial w} &= \frac{\partial \left( 1/2\sum_{n=1}^{N}\{y(x_n,w)-t_n\}^2+\dfrac \lambda2||w||^2 \right)}{\partial w} \\ &= \frac{\partial \left( 1/2(y({\bf x},{\bf w} )-{\bf t})^T(y({\bf x},{\bf w} )-{\bf t})+\dfrac \lambda2{\bf w}^T{\bf w} \right)}{\partial {\bf w}} \\ &= \frac{\partial \left( 1/2(y({\bf x},{\bf w} )-{\bf t})^T(y({\bf x},{\bf w} )-{\bf t})\right)}{\partial {\bf w}}+\frac{\partial \left( \dfrac \lambda2{\bf w}^T{\bf w} \right)}{\partial {\bf w}} \\ &= \left( \frac{\partial y({\bf x},{\bf w} )}{\partial {\bf w} } \right)^T \left( y({\bf x},{\bf w} )-{\bf t} \right) + \lambda{ \bf w }\end{align}$$


Numerical error while implementing a periodic pulse signal

Your "numeric trick" is not working. After running above code, have a look at :

plt.figure(); plt.plot(comb); plt.xlim((324,336)); plt.ylim((-0.1, 1.1))

Your "Dirac" train has peaks with width 2, hence you have the double-height amplitude pulses.

You can fix this by e.g. doing

comb2 = np.zeros(npts)
step = 1e-9*fs  # distance between peaks in terms of samples (floating point)
peakPos = step*np.arange(int(npts / step))  # integer position of the peaks
comb2[np.round(peakPos).astype(int)] = 1  # set the peaks
plt.figure()
plt.plot(comb)
plt.plot(comb2); plt.ylim((-0.1, 1.1)); plt.xlim((324, 336))


signal sampling for 2nd time of sampled signal

thirI believe it is a nonsensical question, or at least a tricky, meant  to detect students who have  a shaky knowledge, and whose anwser deserves a lot of care.

First trick: if the $600$ Hz signal is sampled at  $1000$ Hz, it will endure aliasing (a proper sampling frequency would be above  $2 \times 600$ Hz). If the sampling process is perfect, it will appear as a sine at the mirror frequency $2*\frac{1000}{2}-600 = 400$ Hz.

The  signal (that is assume real)) is now  in the $[0,500] $ Hz band. An ideal-low pass cutting at $500$ Hz would thus be an all-pass, second trick, and would be utterly  hard to design, numerically, in that case (third trick).

So in theory, the  output signal would be the same (aliased) signal as the one from the first stage, fourth trick in a way.


Detecting characters in high contrast/oversaturated regions

Use your camera's RAW mode, if it has one, and underexpose. If these concepts are not familiar to you, head over to https://photo.stackexchange.com/

Is the first picture a close-up of the second one, because the letter K on the tag is quite legible? If not, you may also want to use a higher resolution camera, and/or not shoot from such a distance.

If you think about it, you are asking us to robustly encode roughly four bits (on letter of the English alphabet), judging by the first picture. You could simply map the binary representation of the letter to five boxes; four in a square and one on the side, so you don't get any confusion due to symmetry. (If you laid them side by side you wouldn't know which side is which.) If you need to encode multiple letters, you could try using QR or a barcode.

Your first order of business has to be getting the exposure under control. I'd investigate the possibility of using aftermarket firmware if the manufacturer's doesn't give you enough control.


See pages 5 and 6 and the plots on the following pages specific to number of samples per symbol in this very helpful reference by Ken Gentile on designing RRC pulse shape filters:

http://www.analog.com/media/en/technical-documentation/application-notes/AN-922.pdf

An example I have previously done shows the consideration of filter length (how many symbols does the filter span) and number of samples per symbol. The top plot shows a 100 tap RRC filter with only 2 samples per symbol and the lower plot it is the same number of taps but 10 samples per symbol.

Observe the following key points that illustrate some of the important considerations involved: With 2 samples per symbol the filter impulse response extends further in time. This results in a greater stop band attenuation for the same filter roll-off factor (Alpha) and complexity (number of taps). The difference is particularly dominant close to the band edge in the area of 700 Hz in the plots.  In the figures shown we see that we get approximately 70 dB of attenuation for 2 samples per symbol but only approximately 35 dB of attenuation for 10 samples per symbol. However the overall frequency of the digital spectrum is most constrained for 2 samples per symbol; this is resolved with subsequent up-sampling or consideration elsewhere in our overall system design. As long as those considerations are properly addressed, according to Nyquist sampling theorem we will achieve the full integrity of our RRC filter design (with only 2 samples per symbol!).
In the second case, because of the higher number of samples per symbol, with the same number of filter taps (complexity) the impulse response of the filter does not extend as far in time, and is therefore truncated earlier. This degrades overall performance as visible here with the inferior filter rejection. However the greater overall digital frequency spectrum is wider.
An efficient design would choose to do two samples per symbol for such a "shaping" filter and then upsample to the bandwidth needed to be compatible with the analog design or other bandwidth considerations, since the upsampling filter design would be relatively simpler. Also to be considered for many applications is a polyphase RRC upsampling filter structure. In this case we would choose an even number of taps and incorporate the upsampling as described further in this response: How to implement Polyphase filter?

Keep in mind, as explained further in this post: Why root raised cosine filter can eliminate intersymbol interference (ISI) ? the entire motivation to do RRC filtering is spectral efficiency: hence regulatory requirements (spectral mask) will often dictate the rejection ultimately needed in the RRC filter (as long as subsequent analog stages, in particular the power amplifier, maintain sufficient linearity to not undo the shaping we so carefully applied...resulting in "spectral regrowth!).


Determine causality and stability from given filter structure

It looks like a Homework problem so I will only provide small hints:

Check from the block diagram implementation of the system whether the current output sample $y[n]$ depends on only current input sample $x[n]$ and past input & output samples $x[n-n_0]$ and $y[n-m_0]$. If so, the system is causal, else it will be non-causal.

You can use the information from the previous point (causal or non-causal) along with the pole location of $H(z)$ to check if the system is BIBO stable or not.

Any Causal System with poles inside unit-circle are BIBO stable.
You can refer to this Pole-Zero Analysis for better understanding.


Calculating the inverse filter for the (exponential) sine sweep Method

Assuming that your Exponential Sweep Sine was generated using the formula:

$$x(t)=\sin\left(\frac{2\pi f_1 T}{R}\left(e^{\frac{t R}{T}} -1\right) \right)$$

where:

$f_1, f_2$ - Initial and final frequency of the sweep

$T$ - Duration of the sweep

$R = \ln\left(\frac{f_2}{f_1} \right)$ - Sweep rate

Then the inverse filter is calculated by scaling the amplitude of time reversed $x(t)$ by:

$k = e^{\frac{tR}{T}}$

Which will result in an exponentially decaying sweep:

$f(t) = \frac{x_{inv}(t)}{k}$

Example in Python:

#!/usr/bin/env python

from __future__ import division
import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt

def dbfft(x, fs, win=None):
    N = len(x)  # Length of input sequence

    if win is None:
        win = np.ones(x.shape)
    if len(x) != len(win):
            raise ValueError('Signal and window must be of the same length')
    x = x * win

    # Calculate real FFT and frequency vector
    sp = np.fft.rfft(x)
    freq = np.arange((N / 2) + 1) / (float(N) / fs)

    # Scale the magnitude of FFT by window and factor of 2,
    # because we are using half of FFT spectrum.
    s_mag = np.abs(sp) * 2 / np.sum(win)

    # Convert to dBFS
    ref = s_mag.max()
    s_dbfs = 20 * np.log10(s_mag/ref)

    return freq, s_dbfs

if __name__ == "__main__":
    # Sweep Parameters
    f1 = 10
    f2 = 100
    T = 3
    fs = 1000
    t = np.arange(0, T*fs)/fs
    R = np.log(f2/f1)

    # ESS generation
    x = np.sin((2*np.pi*f1*T/R)*(np.exp(t*R/T)-1))
    # Inverse filter
    k = np.exp(t*R/T)
    f = x[::-1]/k
    # Impulse response
    ir = sig.fftconvolve(x, f, mode='same')

    # Get spectra of all signals
    freq, Xdb = dbfft(x, fs)
    freq, Fdb = dbfft(f, fs)
    freq, IRdb = dbfft(ir, fs)

    plt.figure()
    plt.subplot(3,1,1)
    plt.grid()
    plt.plot(t, x)
    plt.title('ESS')
    plt.subplot(3,1,2)
    plt.grid()
    plt.plot(t, f)
    plt.title('Inverse filter')
    plt.subplot(3,1,3)
    plt.grid()
    plt.plot(t, ir)
    plt.title('Impulse response')

    plt.figure()
    plt.grid()
    plt.semilogx(freq, Xdb, label='ESS')
    plt.semilogx(freq, Fdb, label='Inverse filter')
    plt.semilogx(freq, IRdb, label='IR')
    plt.title('Spectrum')
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('Amplitude [dBFS]')
    plt.legend()

    plt.show()

And output:

For Robert, here is the spectrum plot for case without amplitude modulation of the inverse filter:

Relevant literature:

  Q. Meng - Impulse Response Measurement with Sine Sweeps and Amplitude
  Modulation
  Schemes

  A. Novak - Nonlinear System Identification Using Exponential
  Swept-Sine
  Signal

  K. Vetter - ExpoChirpToolbox - a Pure Data implementation of ESS
  impulse response
  measurement


Why bit rate is same in a lower frequency spectrum and higher spectrum?

If you are a bass and sing A 220 Hz and B 247 Hz at 120 beats per minute, you can sing data at a certain rate.

If you are a soprano and sing two quarter notes A 880 Hz and B 988 Hz at the same BPM, your data rate isn't any higher.

If you sing more notes (higher bandwidth) you could communicate a more complex score (e.g. a higher data rate).  So bandwidth (range of notes) matters, not being a tenor or alto-soprano (spectral frequency).

Same at RF frequencies.


Applying a Power on a Sum of Gaussians to Enhance Resolution

The binomial theorem is not necessarily involved: the top waveform is simply pointwise raised to the 5-th power, as others have noted, and there are no missing peaks that should be seen. To illustrate, consider the following crude hand-drawn figure:

I drew this in my simulation software, so it was automatically "digitized" as I drew it. Obviously, I am poor at freehand drawing with a computer mouse, but I was only trying to make two peaks, with no underlying peak shape in mind.

Then the second figure shows what happens if the digitized values from the first figure are raised to the 5-th power and plotted:

Still ugly, but better resolved.

Now taking into account a more recent answer by @MBaz (https://dsp.stackexchange.com/a/60750/41790), there is still one loose end. To see this, consider two Gaussians, $G_1(x)$ and $G_2(x)$, and adding them. Their sum function, i.e., $G_1(x)$ + $G_2(x)$, does not look like either function. Of course, we know the sum function's two constituent summand functions. But the sum function does not show them: it only shows their sum.

Now square the sum function. Then the pointwise squaring works as expected. Alternatively, if we use the expansion, then we get three constituent functions: $G_1^2(x)$, $G_2^2(x)$ and $2G_1(x)G_2(x)$. So the squared sum function is just the sum of these three constituent functions and none of the three are 'visible' to us: we only see the sum of the three constituent functions. So this is why we can never see the 'extra' peak functions. Indeed, we cannot see any of the constituent functions. We can only see the their sum and should not have expected that sum to display its constituent functions.


ARMA Filter Output Stationary and set up?

Hi: The use of the term "stationary" is vague because there are a few different kinds. But usually, the term "stationary" is meant to mean wide-sense stationary. But, since you've also used the term wide sense stationary in your question, I'm thinking that you're use of the term stationary might just be referring to the output ( by output, I'm assuming you mean the estimated residual noise term )  having a constant (zero) mean and constant variance.
Wide sense stationary refers to constant mean, constant variance and $corr(X_{t1}, X_{t2})$ being only a function of, the time difference, $t2 - t1$.
But, the general process in ARMA modelling ( in statistical time series that is. dsp could be different), is to take the estimated residuals and check that they have a  constant mean of zero, constant variance and that they are uncorrelated. This is because the ARMA framework assumes that the error term is white noise so the estimated residuals should be uncorrelated and have zero mean and constant variance.
So, the test that is done after the ARMA model is identified and estimated
is testing that the estimated residuals are white noise. But, white noise is a specific case of stationarity so one could argue that you're kind of testing for wide-sense stationarity ( except that corr should be zero so not even a function of the time difference ).
Note that there is another term called "strict sense stationarity" which is a stronger form of stationarity where one requires that the joint distribution of  $X_{1}, \ldots, X_{n}$ is the same as the joint distribution of  $X_{t+1}, \ldots, X_{t+n}$.