Combine FFT's of shorter length than sample data to get spectrum of all data

Yes you can do it, but you will loose spectral resolution. The method you mention is a variant of the Welch's spectogram estimation algorithm which relies on the block based estimation of power spectrum of a long sequence in terms of shorter blocks within it.

Given that you have a data of $100$k in length which was sampled at $F_s = 100$k, then a full length DFT after rectangular windowing would provide you a spectral resolution of about $1$ Hz to $3$ Hz, where as a window of $4096$ samples would yield about $25$ Hz to $75$ Hz of spectral resolution.

If this resolution is enough for your purposes, then you can (and indeed you should better) perform the block based averaged spectral analysis in order to trade-off between the reduced spectral resolution and an increased estimation accuracy (decreased estimation variance)


In this particular case, $\nabla I$ is the gradient because the image $I$ is a scalar field. The summation over $\Omega$ would make more sense if $I$ was a vector field. Since the image pixels are samples over fixed locations in space, the $\Omega$ can be dropped.

To understand why this is used as a fitness function, you need to understand Total Variation Minimisation. The assumption there is that a signal has some broad structure that is disturbed by noise. Noise is seen here as large deviations around a central value. But this view would include edges which do not belong to the "noise" part of an image but to what gives it "structure". So, what you are trying to do with total variation regularisation is to estimate the total length of the signal "curve" which largely depends on its structure and then try to minimise the differences between the signal values in areas that it does not "cost" much to do so. If the denoising algorithm tried to suppress an image edge, that would mean a huge change in the length of the curve. But if the denoising algorithm tried to suppress some spurious spike (as low pass filtering would do), then that would not present a large error.

What the authors are trying to do is to use the outputs of three denoising algorithms to descend to the lowest error of this fitness function as quickly as possible by selecting the "best" areas from each denoised image and blend them together with low pass filtering.

Hope this helps.


Can the overlap-add method be used with varying cut-off frequency filter?

The problem here is not the application of the FFT for filtering but the fact that you want to model a linear time-variant system using methods for linear time-invariant (LTI) systems.

If you used the linear convolution instead of the circular convolution + overlap-add you would still face the problem that you have to divide the audio signal into blocks of data during which the impulse response is constant.

One solution is to make the data blocks sufficiently short so that the impulse response can be considered quasi-invariant during one data block. For example, if your audio sample rate is 44.100 Hz and your FFT length is 1024 every impulse response change needs to be approximated by 43 discrete steps per second.


Relationship between interpolation accuracy, impulse response and frequency response

any interpolation polynomial that goes through the two adjacent samples to the left and right of the area being interpolated will result in an effective impulse response that is a piecewise-polynomial function going through zero at all integers except zero (like the sinc function does).

e.g. linear interpolation:

third-order lagrange:

third-order hermite:

note all of these effective impulse responses go through the $(x,y)$ points of $(0,1)$ and $(n,0)$ for any other integer $n \ne 0$, just like the sinc function does.

so you can represent that interpolation function as something times a sinc function and the "something" can be obtained by dividing that interpolation function, point-by-point, with the sinc function.  the only division by zero is when the numerator is also zero and that can be, what we like to call in the hand-waving department, a "removeable singularity".  (maybe we can find a decent limit at the $\tfrac00$ divisions.)

any interpolation that forces the smooth interpolated function to go through the discrete points can always be described as a windowed-sinc function.  all you have to do is divide to get the effective window.


How to recognize colors from the color chart?

You can house hough lines to detect the color area

First detect all lines
Lines with longest length would be the outer most ones
Pick the top horizontal line and do template match for the squares.
Once you find the squares, traverse downwards toll you reach the other set of squares.
Mark all points where the squares were found
You now have an enclosing area for the color chart.


How do I convert a filter from Spatial to Fourier domain?

For a low-pass filter that is small compared to the image you can put your filter (unscaled) in the center of a square/rectangle the size and shape of your image (all the other pixels are 0), take the FFT and multiply the result  in the Fourier domain.

Since you are using a Gaussian filter you could do something more efficient than that.  The Gaussian filter is _ separable_ so you can actually create a 2D Gaussian filter by doing a 1D Gaussian filter in every row (column) in one dimension followed by a 1D Gaussian filter in every column (row) in the other dimension.

Second, the Fourier transform of a Gaussian is a Gaussian.  So we can work in one dimension at a time, we know what we're going to get, and multiplying a Gaussian in the Fourier domain is also separable since $e^{-(x^2+y^2)}=e^{-x^2}e^{-y^2}$.  The only question is what the standard deviation ($\sigma$) of the transformed Gaussian is going to be.  The derivation is worked out on Mathworld as:
$$
\mathcal{F}[e^{-ax^2}](k)=\sqrt{\frac{\pi}{a}}e^{-\pi^2k^2/a}.
$$


Once you remove the systematic error (or trend) you want to test whether the process is zero mean or not. In hypotheses testing this would assume the following formulation:

$H_1: \quad x[n] = \theta + w[n], \quad n=0, \dots, N-1$,

$H_0: \quad x[n] = w[n], \quad n=0, \dots, N-1$,

where $w[n] \overset{i.i.d}{\sim} N\left ( 0, \sigma^2 \right )$, $\sigma^2$ is an unknown nuisance deterministic parameter and $\theta$ is another deterministic parameter that bares information on determining which hypothesis we choose.

We will develop a generalized Likelihood Ratio Test (GLRT) for this problem. In the GLRT framework the unknown deterministic parameters are estimated using maximum likelihood criterion, i.e.:

$\quad \widehat{\theta}_{ML} = \frac{1}{N}\sum\limits_{n=0}^{N-1} x[n]$, is the sample mean.
$\widehat{\sigma}^2_{ML} = \frac{1}{N} \sum\limits_{n=0}^{N-1} \left ( x[n] - \widehat{\theta}_{ML}  \right )^2$ is the sample variance under $H_1$, and,
$\widehat{\sigma}^2_{ML} = \frac{1}{N} \sum\limits_{n=0}^{N-1}  x^2[n]$, under $H_0$.

After some algebra dividing the two likelihood functions we get the following test:

\begin{equation}
Z = \frac{\widehat{\theta}_{ML}^2}{\hat{\sigma}^2_{ML}|H_0} = \frac{\left (  \frac{1}{N}\sum\limits_{n=0}^{N-1} x[n]\right )^2}{\frac{1}{N} \sum\limits_{n=0}^{N-1}  x^2[n]} \overset{\underset{>}{H_1}}{\underset{H_0}{<}} \gamma,
\end{equation}

where $\gamma$ is determined by the Neyman-Pearson's theorem given a desired probability of false alarm $P_{FA}$.


Since in your system there is no aliasing and you do not filter away important parts of the signal (because LP is having the correct cutoff-frequency), there is no information loss.

When a system changes the signal shape, it does not necessarily mean that information is lost. As long, as the operation is reversible, the same amount of information is contained in $y[l]$ and $x[n]$. However, even though the same information is contained, it can be differently difficult, how to extract the useful information from a given signal, depending on its shape.

However, in reality this ideal low-pass filter does not exist, and it might lose some information.


Formula of frequency modulated sine wave from data points

If your function is $y(x)=\sin(2\pi fx+C)$, and you have data points $\{x_i,y_i\}$, then finding $C$ is easy, if you suppose $x=0$ then $C=\arcsin(y(0))$.

Now, according to your plot your frequency varies with $x$. The first thing you need to do is to know how does your frequency vary with $x$. Supposedly this variation is linear, then you write: $f(x)=ax+b$ and your function model becomes $y(x)=\sin(2\pi f(x)x+C)$.

The problem now is finding (an estimate for) $f(x)$. You have two unknowns to determine $a$ and $b$, you need two data points ($\{x_1,y_1\}, \{x_2,y_2\}$) and you solve a system of two (nonlinear) equations with two unknowns.

$$\begin{cases}
y_1=\sin(2\pi (ax_1+b) x_1+C)\\
y_2=\sin(2\pi (ax_2+b) x_2+C)
\end{cases}$$


Design of FIR Filters with Arbitrary Magnitude and Phase Responses

For the sake of clarity, let me point out that you do not approximate $(1)$ by $(2)$, but the other way around: you approximate a desired frequency response, i.e., your specification, by the filter's frequency response, and this is done by choosing the coefficients in the vector $\mathbf{h}$ such that some error measure (e.g., least squares) is minimized.

Since the frequency response $(1)$ is linear in $\mathbf{h}$, the problem is a linear least squares problem, which can be solved by solving a system of linear equations. This is true for the continuous error measure $(3)$ as well as for the discrete error measure

$$E(\mathbf{h})=\sum_{k=1}^K|H(e^{j\omega_k})-D(e^{j\omega_k})|^2\tag{1}$$

where $K$ is the number of discrete frequency points $\omega_k$. In general we require $K\gg N$, where $N$ is the number of FIR filter coefficients in the vector $\mathbf{h}$.

In the general case it is more convenient to use the discrete error measure $(1)$. There are only a few special cases in which the integral of the continuous error measure (Eq. $(3)$ in your question) can be solved analytically, otherwise you have to resort to numerical approximations anyway. And even if an analytical solution is available, the difference when minimizing the discrete error measure is negligible if the number of frequency points $K$ is sufficiently large.

Minimizing the discrete squared error $(1)$ is equivalent to solving the following overdetermined linear system of equations in a least squares sense:

$$H(e^{j\omega_k})\stackrel{!}{=}D(e^{j\omega_k}),\qquad k=1,2,\ldots,K\tag{2}$$

Formulation $(2)$ is convenient because many software packages such as Matlab have syntactically elegant and computationally efficient methods to solve such systems.

I've implemented a complex FIR filter design based on solving $(2)$ in the (extremely simple) Matlab/Octave function cfirls.m.


Is there any special when all Fourier components have the same angle?

sine waves all have zero crossings at the same phase of $0$ or $\pi$ and have a peak (or valley) with a phase of $\pm\tfrac{\pi}{2}$.

cosine waves all have zero crossings at the same phase of $\pm\tfrac{\pi}{2}$ and have a peak (or valley) with a phase of $0$ or $\pi$.

so if it's a $\sin(\cdot)$ series (with no cosine terms), expect a big zero-crossing when all of the phases line up at $0$ or $\pi$.

or expect a $\cos(\cdot)$ series (with no sine terms) to have a major zero-crossing transition to occur when the phases are $\pm\tfrac{\pi}{2}$.


Can two nonzero signals $x[n]$ and $y[n]$ give a zero convolution

Yes, for example let

$$x[k]=1$$

for all $k$ and

$$y[k] = \begin{cases}1 & k=0\\-1 & k=1\\0 & otherwise \end{cases}$$

It is easy to see that in case of a convolution, the result will be zero for all values of $n$.


Optimising an input signal for a system to minimise distortion of output signal

System is non linear so all bets are off.

Jokes aside though, wouldn't you just characterize the system, then generate an inverse system that the signal passes through first such that the effects of the original system are cancelled out?

Edit: another potential solution is to 'linearize' the system by finding the characteristics of it, then scaling your input signal such that the system behaves in a linear way. Then follow the inverse idea that I gave in my original comment.


What distribution is the easiest to compress?

It's very easy to make fundamental mistakes on the topic of information theory, unless you have some PhD or alike in the field. That being said as far as I remember the following could roughly be stated about your question.

Considering a discrete probability and source model of iid, the most efficient compression of a source happens when the source has minimum Entropy (information) in it; i.e, its average codeword length (in bits) per symbol is minimum. Such a variable length code would have the highest compression and the most efficient coding to represent the source in a lossless manner.

Such a case would happen when one of the symbols of the source alphabet has a probability close to 1, and all other probabilities are close to zero. Such a distribution is said to be highly skewed. But I'm not sure if a unique minimum can be shown to exist unlike the case of maximum Entropy (minimum compression) which happens when all the symbols have equal probability.

On the other hand, you are talking about lossy Quantizers and probability distributions. PDF optimized lossy quantizer design is about finding the best input decision and output reconstruction levels for a given soure PDF such that the the expected quantization error (MMSQE) per given source distribution a minimum. And that's the Max-Lloyd quantizer for the scalar case. Rate-Distortion theory shall be invoked for various considerations here.

Merging these two concepts you would say that, for a given source power (or dynamic range) the one which required minimum number of quantization levels while still providing least possible quantization error (distortion) would be given by a source whose PDF is dense at a certain region. In that case you would put a few fine quantizer levels at those dense signal positions and then use larger step size levels outside of that region, which would still provide minimum expected distortion.


A few notes regarding your approach and detailed questions:

First, it is very common in audio analysis to split the signal to be analyzed into short overlapping windows, because audio signals are not stationary (their characteristics change over time), so they need to be processed on short segments over which they can be considered stationary (many analyses assume stationarity). Typical segment size is in the 5 .. 20 ms range. You'll be left with an odd-size segment at the end, just pad it with zeros and that's it.
In speech recognition the most widely features are MFCC. AR coefficients are quite common in older literature too. Zero crossing rate or energy are unlikely to discriminate phonemes in a robust way.

Now, regarding your project...

First, it is not clear if the goal of your project is to arrive at a working system even if you have to reuse existing code for that, or to do it yourself to learn on the road...

If you want to reach the goal without learning anything, there are ready to use speech recognition packages like Sphinx or HTK - the catch is that to understand these tools one actually has to know a lot about speech recognition tech. The setup time (training plans, model topology, grammar) and the concepts can be daunting.

An intermediate level would be to use third party code for signal analysis and feature extraction, and then focus on the recognition itself. You'll find C/C++ MFCC or AR modeling code in audio feature extraction packages like YAAFE or Aubio. You can also use the feature extraction command line utility of HTK (HCopy with the right config file for extraction + HList to get a clear text file).

When designing a speech recognition system, you have to ask yourself these questions:

What is the vocabulary size (from discriminating 2 words to recognizing natural language)?
What is the speaker variability (will I have to recognize from the same speaker again and again, or will I have to deal with many different speakers)?
What kind of training data will I have, in particular, is it possible to collect labelled speech samples from the same voice as the one you'll have to recognize?
What are the conditions (clear recordings... or voice through a glitchy mobile phone with background street noises)?

You have only answered the first question. If the answer to the other ones are these ideal conditions: the system will only deal with one speaker + will have a corpus of a few recordings of yes/no uttered by this same speaker + there will be no background noise ; the problem is not that hard. And these ideals conditions are not that unrealistic - there are some voice control PC applications which require this kind of prior training, the user needs to pronounce a few occurrences of each command word.

In these ideal conditions, a rustic approach (early 80s tech), but which will work, is DTW. Extract from each utterance in the training corpus (which will be labelled either "yes" or "no") a sequence of feature vectors, and compute the DTW distance between the utterance you need to recognize and those in the training corpus. Use this distance metric and a "k-nearest neighbour" scheme to decide on the label. DTW is archaic compared to a proper acoustic model + word model + language model graph, but for a very simple vocabulary (just a handful of words) it does the job.

Now, if you want something speaker independent (which will work for many speakers, including those for which we don't have a template), you'll have three plans. One is to stay with the same approach and record utterances from a very diverse set of speakers in terms of gender, accent, hoping that it will cover enough ground. But DTW scales very poorly! Another is to use an adaptation scheme - find the transform of your feature vectors that minimize the sum of DTW distance (or the minimum DTW distance) between your utterance and the database (paradoxically, I'm sure you'll find more literature about adaptation in the context of classic HMM-based recognition instead of DTW). The last approach is to forego DTW and do the real thing (phone-level GMM + word + language model).

I won't be as negative as the other members. I used to give a very similar project (digit recognition) to students which had no prior experience with audio (but a semester of signal processing and prior experience with Matlab + some background literature before starting the project), and they completed it in Matlab in two or three 3 hours lab sessions (roughly one for MFCC computation, one for DTW, one for the whole corpus management and k-NN + evaluation routines) with about 100 lines of code in total - so this is not an unimaginable engineering feat requiring a full team of PhDs...


Audio Watermarking in MATLAB - How to choose frames to embed in?

Loudness is defined as human perception of volume. In general, there is a very good correspondence between the two. Therefore, you can often use one as a proxy for the other.

In this case, you can easily measure volume: the portions of the audio signal with the highest signal intensity (or volume) correspond to the parts with the greatest absolute value (or, if you prefer, square value). With a few assumptions, you can say that the loud sections are also the sections with high volume.


In my opinion, the most straightforward way to solve this problem is to consider the time domain signal as it passes through the signal chain. A unit step at the input when sampled and put through a ZOH becomes just another continuous-time unit step. Consequently, the Laplace transform of the continuous-time signal just before the sampler at the output is given by

$$C(s)=\frac{e^{-s}}{s(s+a)}=\frac{e^{-s}}{a}\left[\frac{1}{s}-\frac{1}{s+a}\right]\tag{1}$$

which corresponds to the time domain signal

$$c(t)=\frac{1}{a}\left[1-e^{-a(t-1)}\right]u(t-1)\tag{2}$$

Sampling $c(t)$ at $t=kT$ gives

$$c(kT)=\frac{1}{a}\left[1-e^{-a(kT-1)}\right]u(kT-1)\tag{3}$$


Properties of a custom wavelet family for continuous wavelet transform

It sounds like you don't need a wavelet. A wavelet's translation and dilation properties allow you to get both time and frequency resolution.

A matched filter takes an inner product of a signal with another signal, allowing you to detect a known signal in noise (by correlating the signal you have with the known signal). You may want to use this idea instead.

Another idea is matching pursuit type algorithms -- if you have an overcomplete set of signals, and you want to find an approximation of a signal with respect to these overcomplete set of signals, you can do this through a matching pursuit algorithm. Basically, iteratively form a residual by finding the largest inner product signal in the set and then subtracting a scaled verison of the signal.


Calculating phase response of filter with latency

Thanks to everyone who provided helpful comments. As some suggested shifting the signal might be the way to go, but this results in the blue line where the green line is the phase I would expect.

It turns out to correct the phase after taking the DFT of the impulse response of the filter when the filter has a delay each bin needs the following angle adding to it:

$$
\frac{ 2\pi i L} {N}
$$

Where $L$ is the latency in samples, $N$ is the FFT size and $i$ is the $i$th FFT bin. So I was missing this as step 5.

For anyone else who stumbles upon this, the phase will end up being neither wrapped or unwrapped. So you will need to take the remainder of the final phase modulus $\pi$. This will give you the normal wrapped phase (as a correctly quadranted arctan would produce by just taking the phase of the complex valued FFT bin) and then you can unwrap as normal.


How to proceed for BOW creation From sift features

There is now support for the bag-of-words model in the Computer Vision System Toolbox for MATLAB.


What is the value of the standard deviation $\sigma$?

randn gives you Gaussian noise with std. dev. of 1.

(x/255)*(randn(size(image))) gives you a noisy image of std. dev. (x/255)


Here's an attempt to do what you need in scilab.

The top plot shows some data that I synthesized.  The second plot shows the auto-correlation of the raw data (zoomed in around the central peak of the auto-correlation).  The red circles show the peaks found using this find_peaks function.

The final plot shows the difference between all the peak locations. This will be an estimate of the period.  Because you're not really guaranteed that the underlying period will be an integer number of samples, you should probably find the mean of these values.

In this case, the "true" period is 1/f0 = 11.191996, and taking diffs = diff(peaks); and then mean(diffs(10:173)) yields 11.195122.

Code below.

N = 1000;
f0 = 0.0893495634;
phi = rand(1,1,'uniform')*2*%pi;
sigma = 0.5;

x = sin(2*%pi*[0:N-1]*f0 + phi) + sigma*rand(1,N,'normal');
XC = xcorr(x);

clf
subplot(311)
plot(x);

subplot(312)
plot(XC);

peaks=peak_detect(XC,0);

plot(peaks,XC(peaks),'ro')

a = get('current_axes');
a.data_bounds=[950 1050 -500 800];

subplot(313)
plot(diff(peaks));


You will need to use a continues wavelet transform (cwt) where you can specify the resolution of the scales (as done in the image you posted mathworks.fr/help/releases/R2013b/wavelet/ref/cwt.gif). The DWT only has $J$ scales by definition for $2^J$ signals.


These are just different ways to implement FIR filters. In theory these structures are all equivalent in the sense that they compute exactly the same output signal for a given input signal (under ideal circumstances, i.e., infinite precision), but in practice they behave differently. The differences are mainly the following:

their behavior when implemented with finite word-length
computational complexity
memory requirements

I'm sure that if you continue reading the book, these structures will be explained in more detail, and their relative merits will be discussed.


What are physical scenarious which can be modeled as multiplicative noise?

In communications, a very important multiplicative model is a flat fading channel where the received signal can be modeled as

$$y(t)=h(t)\cdot s(t)+n(t)\tag{1}$$

where $h(t)$ is the time-varying attenuation caused by the channel, and $n(t)$ is additive noise.


Encounter 0 when calculating log power spectrum

Typically, one's data has some sort of noise floor (-96 dB, etc.).  So one common way to deal with the FFT bins that are zero (or tiny) is to replace anything below the noise floor with the noise floor level, since any value below that value is most likely not useful data.  Doing this before taking the log() function may provide some computational efficiency.


How to derive a rate (e.g. water dripping) from an audio file

Beat detection is not what you are looking for.

Beat detection aims at looking for a periodic, sometimes implicit pulse in music signals. The events you are interested here are not periodic and every drip counts. What you are looking for is instead "Onset detection", to recognize each percussive event in the audio signal. There are ready to use software solution like aubio or sonic annotator to extract the timestamps at which percussive events occur in an audio file. From there, it is only a matter of running a sliding window over your output, count events within the window and divide by the duration of the window.


If a filter has a non-linear phase lag will the output signal be distorted?

I think your question comes from several misunderstandings. The fact that the phase lag of a system becomes more negative for large frequencies does not mean that there's more distortion of larger frequencies. Neither does it mean that high frequency signal components experience more delay when passing through the filter.

Imagine an ideal system that simply delays the input signal while keeping its amplitude unchanged. So for an input signal $x(t)=\sin(\omega t)$ the output signal is given by

$$y(t)=\sin(\omega(t-\tau))$$

where $\tau$ is the delay. The output signal $y(t)$ can be written as

$$y(t)=\sin(\omega t-\omega\tau)=\sin(\omega t+\phi(\omega))$$

where $\phi(\omega)=-\omega\tau$ is the phase response of the system. Note that the phase becomes more negative for larger frequencies, and this is exactly what is needed for a system that does not distort the input signal.

What causes distortion are deviations of the system's magnitude response from the value $1$, and deviations of the system's phase response from the ideal linear phase shown above.

And, as mentioned in the comments, infinite impulse response (IIR) filters like a Butterworth filter are almost always implemented in the time domain. You implement a given difference equation, which is equivalent to implementing a system with a given impulse response and a given frequency response (the latter two are Fourier transformation pairs), and the phase response is part of the system's frequency response.


Problem with calculating recursive equation

Sorry for answering like this but I cannot comment since my reputation is not sufficient, but yet I want to help.

It is a filter and you can do it way easier, and more quickly than just building an algorithm(in Matlab ofcourse). It will take relatively a long post to cover all the concepts in order to make a sufficient understanding about filter designs.

However, I can help you about which topics you should cover.

Please take a look at FIR/IIR filters (Your recursive equation is basically an IIR filter) and Z-Transform topics.

If you need a quicker way to learn about IIR filter this lecture will give you some information.


Is an image a random variable or a random process?

While there are some obvious exception to it ("static" pattern on a television screen, the "dark frame" noise pattern of a camera), images are rarely generated by random processes. Declaring that an image is drawn from such or such distribution or generated by such or such random process is just a post-hoc modeling decision, and there is no "ground truth" to validate or invalidate this choice; other than the performance of the machine vision / image enhancement etc. method derived from the modeling decision.

So you can view a $m \times n$ image as a single random matrix (I assume this is what you mean by your first alternative - the image as a whole being considered as a single multidimensional random variable); or you can view it as a random field (a collection of random variables indexed by $|[1, M]| \times |[1, N]|$). I've encountered the random field view more often than the random matrix view.

When using a random-field approach, you can view each pixel as i.i.d ; or you can introduce dependencies between pixel values through the use of a Markov Random Field model. These are not the only option - you could very well consider a two-layer model where one first random process assigns a region index to each pixel of the image, and then where the value of each pixel is drawn from a distribution indexed by the region id! No approach is "better" than another. The more complex the model, the more "plausible" the images it will generate, but the more intractable the computations might become. It often helps, when using this kind of statistical approaches, to draw a few sample images from the chosen distribution/random process, and look at them to have a good grasp of what kind of assumptions you have built into your model.


Find $X(j\omega)$ after sampling of $2\cos(2000\pi t)+\sin(5000\pi t)$ at 5 kHz sampling rate

Look at your signal and figure out what sampling rate you need to avoid aliasing. In this case, you are sampling at twice the maximum frequency component so you should not expect aliasing, but you work leads to aliasing components showing up at plus/minus $\frac{3\pi}{5}$.

Remember that the DTFT is periodic with period $2\pi$, and try re-working the problem.


Implementing an offline guitar monophonic pitch detector

Each of the "standard methods" you listed has a failure mode (or modes) and an estimation error (which may vary with the guitar and the note). Characterize the failure modes or each method, and evaluate the estimation errors of the more successful methods, then match them to your samples, and perhaps pick the least bad (or other suitable statistical fit.)

With "the guitar in particular", one thing to watch for is that some notes on some guitars can produce a slightly (but measurably) inharmonic overtone series.  If so, the waveform will never be perfectly periodic.  This can cause successive periods of the fundamental mode waveform to never exactly match (which can throw off both autocorrelation and AMDF lag estimation methods, etc.).  Thus, one may need a psychoacoustic model to determine what might likely be the perceived pitch to a human (assuming that's what you care about).  Also, the ratios between overtone magnitudes will usually change (evolve) over time, not just decay linearly as per an ADSR envelope model, which might also affect periodicity and your choice of method.


Determine whether a signal is increasing or decreasing over time

You may use a bandlimited differentiator (i.e. the cascading of a low-pass filter to remove noise, and a differentiator to get the value of the slope).


How to express STFT and ISTFT as a 1d convolution and 1d deconvolution in tensorflow/keras

I found an implementation of STFT based on conv1d in pytorch here:

https://github.com/huyanxin/phasen/blob/master/model/conv_stft.py

edit: Actually, the phasen repository took the STFT code from https://github.com/pseeth/torch-stft
edit2: Asteroid has an alternative implementation of STFT and iSTFT: https://github.com/mpariente/asteroid/


Determining effect of powers of DFT matrix conjugate on input without multiplying

In both cases (odd and even powers of 2) you should get as a result simply a time inversion of your original signal. Because of the implicit periodicity you get from the original vector

$$[x_1,x_2,x_3,\ldots,x_N]$$

the following time-inverted (and periodically continued) vector

$$[x_1,x_N,x_{N-1},\ldots,x_2]$$

Here's an 'odd power' example:

N=8;F=dftmtx(N);U=F/sqrt(N);
x=[1 2 3 4 5 6 7 8]';
y=U'*U'*x;
y =

1.00000 - 0.00000i
8.00000 + 0.00000i
7.00000 + 0.00000i
6.00000 - 0.00000i
5.00000 + 0.00000i
4.00000 + 0.00000i
3.00000 + 0.00000i
2.00000 - 0.00000i


Mathematical justification for zero padding?

Resolution (in the context of spectral analysis) is the ability to distinguish two or more closely spaced sinusoidal components in a spectrum. If you can't resolve them, this means that you only see one maximum instead of two or more. Resolution is determined by the window length, i.e., by the number of time domain samples (given a certain sampling rate). Windowing will smear out or broaden any narrow band components in the spectrum. Remember that windowing, i.e., multiplication in the time domain, corresponds to convolution in the frequency domain, and this convolution with the spectrum of a standard window is a sort of averaging operation in the frequency domain, hence broadening any (frequency domain) impulses.

Resolution can only be improved by adding more data points, i.e., more information, to your analysis window. Zero-padding obviously doesn't add any information. What zero-padding does is increase the sampling density of the discretet-time Fourier transform (DTFT) of the windowed signal. For a length $N$ signal, the DTFT is given by

$$\tilde{X}(\omega)=\sum_{n=0}^{N-1}x[n]e^{-jn\omega}\tag{1}$$

The discrete Fourier transform is

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

Comparing $(1)$ and $(2)$, we see that the DFT is just a sampled version of the DTFT:

$$X[k]=\tilde{X}\left(\frac{2\pi k}{N}\right)\tag{3}$$

Zero-padding increases the number of equidistant samples taken from the DTFT $\tilde{X}(\omega)$.

A detailed analysis of the fact that zero-padding corresponds to interpolation in the frequency domain is given in this answer.


Determine Parameters of gaussian distribution to fit a certain Image-Filter

The filter you referenced is known as a Binomial filter. It is an approximation to a Gaussian, but for smaller filters it's a very crude one. The design is oriented more toward efficiency than accuracy. Rather than sampling a Gaussian directly, the idea behind the approximation is based on the Central Limit Theorem. In this case, it means that a small moving average filter $[ 1 \space  1 ]$ convolved with itself over and over will become more and more Gaussian. It also happens to yield the rows of Pascal's Triangle. The coefficients for larger filters are given by the following formula:

$$
{N \choose k} = \frac{N!}{(N - n)!\space N!}.
$$

For example, for $N = 5$, taking the outer product of the corresponding row of Pascal's triangle with itself yields

$$
\begin{bmatrix} 1 \\ 4 \\ 6 \\ 4 \\ 1 \end{bmatrix}\begin{bmatrix} 1 & 4 & 6 & 4 & 1 \end{bmatrix} =
\begin{bmatrix}
1 &   4 &   6  &  4  &  1 \\
4 &  16 &  24 &  16  &  4 \\
6 &  24 &  36 &  24  &  6 \\
4 &  16 &  24 &  16  &  4 \\
1 &   4 &   6 &   4  &  1
\end{bmatrix}.
$$

The paper I referenced in the first link has the following diagram which shows how the quality of the approximation increases with increasing filter order:


The Distribution of Filtered Gaussian White Noise

If you filter a Gaussian random process with an LTI system, the output will also be Gaussian. You can make intuitive sense of this by considering that a linear combination (which is what filtering does) of jointly Gaussian random variables is a Gaussian random variable.

You can find an in-depth treatment of filtering random processes in this MIT OpenCourseWare document (section 7.4).


What is the point of smoothing an FFT or spectral density plot, and how does that affect the noise floor?

This question is specific to smoothing samples in the frequency domain (given by FFT and spectral density) and asking about the impact to the resulting noise floor in the same domain.
The answer depends on the characteristics of the noise and any time domain windowing that is applied. For white noise, with no windowing beyond the rectangular window, each bin of the DFT is independent of the next and with a moving average across the bins the standard deviation of the noise is reduced by $\sqrt{M}$ where $M$ is the number of samples in the average. The standard deviation of a single tone however would be reduced by $M$, and thus the ratio of the two would go down $\sqrt{M}$.  We may think we are reducing the noise by smoothing the spectrum, but we are reducing the SNR since the signal components of tones (that occupy single bins) will drop more!!
This makes sense intuitively as after the smoothing through a moving average each frequency bin now includes the noise of the adjacent samples in the average.
This is clear from the additive property of independent identically distributed random variables:
The variance $\sigma^2_M$ of $\sum_{k=0}^{M-1}X_k$ is $\sigma^2_M = M\sigma^2_k$ where $\sigma^2_k$ is the variance of each $X_k$.
While for correlated variables (as we would have with tones) the relationship is $M^2\sigma_k^2$. (The signal increases in magnitude at rate M, while the noise increases in power at rate M).
A very simple example may make this clearer: consider the bell curve distribution of white noise process with a non-zero mean: if you added two independent samples from this process, the mean would go up by a factor of $2$ but the standard deviation would only increase by $\sqrt{2}$, but if the noise samples were dependent or specifically the same, then both mean and standard deviation would both increase by $2$ (a simple scaling of the random process).
So if we are trying to access/characterize a noise floor in the presence of a strong signal, this would be recommended to smooth the variability of that noise floor (equivalent to adjusting the Video Bandwidth rather than Resolution Bandwidth on a spectrum analyzer which basically serves to average the noise floor rather than reduce it). But if we are trying to detect a weak signal in the presence of noise this would not be recommended. Continuing the spectrum analyzer analogy, in which case we would reduce the Resolution Bandwidth (RBW) which will result in reducing the noise relative to our signal... for the DFT this means increasing the number of samples, just as in the spectrum analyzer the sweep rate must reduce when we reduce RBW -- time must increase!).
Windowing in time will reduce this change in SNR since the window will impose correlation on adjacent samples (with its own impact on SNR due to that). Since windowing is a product in the time domain which is a convolution in frequency, we can see how these two effects (windowing in time, moving average in frequency) are one and the same when the time domain window is the (aliased) Sinc function, aka the Dirichlet Kernel.


If you're trying to identify birds, you should look into the Mel-frequency Cepstral Coefficients (MFCCs). They are commonly used in speech recognition, and there's actually a fair amount of academic papers on bird song recognition using the MFCCs and various other features.

Short-term energy is definitely useful to extract areas of bird sound, but it probably won't be sufficient to just take an FFT and look for coefficients. For example, imagine the same two  bird sounds, but one is slightly higher pitch, or there is more/different noise present. The FFT coefficients will be significantly different. The MFCCs basically capture the shape of the spectrum, on a scale that emphasizes regions of the spectrum where more detail is needed, and discard some of the unimportant details.


Why complex numbers are represents as a+ib and can't be as (a,b)?

Yes, in signal processing, complex numbers are usually visualized on the complex plane, as you have said.

The reason is that if you put them on a plane, then you are able to measure two important quantities:

1) Magnitude, which is $\sqrt{x^2 + y^2}$

2) Phase angle between your point and the origin, given by $\tan^{-1} \frac{y}{x}$.

If you simply left them as a point, ($x$,$y$), you would not be able to concretize and have a frame work for those quantities.

You may ask, why are those quantities, in turn, important? In signal processing, we are of course dealing with signals, and physically, we are dealing with 'real' signals. However, though a nice trick, an constant oscillation of a quantity in 'real' life, (like a cosine wave), is equivalent to two phasors, rotating around in opposite directions on the complex plane, and adding up together. With this framework, we can see that the phase angles 'cancel' each other out, and that the magnitudes of their resultant give us the magnitude of our 'real' signal.

In fact this is what one of euler's formulas captures. That is:

$$
\cos(2\pi ft) = \frac{e^{j2\pi ft} + e^{-j2\pi ft}}{2}
$$

You can see here how we can easily relate a 'real' world concept, like an oscillating cosine wave, with the 'complex' world of phasors, as they exist and rotate around in the complex plane.

This is one of the corner stones of DSP.


What is the relation between PCM amplitude values and the values in a Spectrogram?

Spectrogram contains magnitude values only, that explains why the values start at 0 instead of some negative value. And scipy seems to calculate it with 24bit accuracy, where $2^{24}\approx 16.7\cdot10^6$.


You need to do it numerically, as a related toy problem $2^a + 3^a = 4$ appears to have no symbolic solution. Bisection method (binary search) is probably the easiest solver to implement:

minexponent = 0;
maxexponent = 2;
targetsum = 10;
precision = 32;
loop precision times {
  exponent = (minexponent + maxexponent)*0.5;
  if (sum(s.^exponent) < targetsum) {
    maxexponent = exponent;
  } else {
    minexponent = exponent;
  }
}
return exponent;

The smallest possible exponent value is 0 for your problem. If you are unsure about the largest possible exponent value, you can test the easy-to-check exponents 1, 2, 4, 8, 16 etc. (as many as needed before the sum goes below the target) before the above code, exponentiating into an auxiliary array by multiplying each element by itself.

Note that if all of the elements are zero, or if all of the elements are one, then you cannot get any other mean than that by exponentiation.

The plot below shows 1000 runs of the algorithm, and the value of exponent at each iteration.


Need help understanding convolution indexing

This question is plagued by bad notation that the OP is unwilling to abandon because, unfortunately, his book uses it. This greatly
impedes his understanding of the concepts.

  Case 1: First convolution equation example:
  $$
h[k] * h^*[-k] = \sum_{n=-\infty}^{\infty} h[n] \ h^*[n+k]
$$

  The way I got the right hand side, was:

  1) Make the $h$'s have dummy variables, so $h[n]$ and $h^*[-n]$.

  2) Flip one of them around, so now I get $h[n]$ and $h^*[n]$

  3) Add the delay $k$ into the one that way flipped, so finally I have $\sum_{n=-\infty}^{\infty} h[n] \ h^*[n+k]$.

${}{}{}{}{}$

No. 3) is not correct. When convolving, you flip one signal and then move it to
the right by $k$ which will give you $h^*[n-k]$, not $h^*[n+k]$.  Thus you end
up with $\sum_n h[n]h^*[n-k]$ as in the leftmost sum in $(2)$ below.

  Great, I have a convolution.  (emphasis added)

The right side of the last equation is not generally called
a convolution; it is called the autocorrelation function of $h$.
If $h$ and $g$ are two finite-energy sequences (meaning $\sum_n |g[n]|^2$
and
$\sum_n |h[n]|^2$ are finite), then the $k$-th term of
their convolution is
$$\sum_{n=-\infty}^\infty h[n]g[k-n] = \sum_{n=-\infty}^\infty h[k-n]g[n].
\tag{1}$$
It is not a good idea to use $h[k]\star g[k]$ to denote the
convolution in general or its $k$-th term in particular
-- far better to write
$h\star g$ in general and $(h\star g)[k]$ in particular -- as the
OP has already been told repeatedly on math.SE. Be that as it may,

  for the case when $g[n] = h^*[-n]$ for all $n$, we can write $(1)$ as
  $$\sum_{n=-\infty}^\infty h[n]h^*[n-k] = \sum_{n=-\infty}^\infty h[k-n]h^*[-n]
= \sum_{m=-\infty}^\infty h[m+k]h^*[m]
\tag{2}$$
  where the last equality was obtained by replacing $n$ by $-m$.
  Notice the difference from the OP's result. For complex-valued
  sequences, there are two different ways of defining the autocorrelation
  function $R_{hh}[\cdot]$ depending on which sequence we choose
  to conjugate, but both ways give the same conjugate symmetry
  property $R_{hh}[-n] = (R_{hh}[n])^*$.

${}{}{}$

  Case3: The problem:

  This is my problem. Case 1 showed the steps for how to do a normal convolution. (emphasis added)
  Case 2 showed how the method of substitution works. However, if I try to solve Case 1 using this substitution method, I get the wrong answer, even though it worked on case 2. Why is that?

  To wit - where am I going wrong here?

Where you are going wrong is that your analysis
of Case 1 is not quite correct
as explained above, while your Case 2 is not a convolution at all,
and is not an autocorrelation or cross-correlation of complex-valued
sequences either, because these require conjugation of one of the
terms.  So you are really mixing up things and then complaining that
the answers are not the same.


Laplace Transform: zeros and corresponding impulse response $h(t)$

For you first question you can use the following
$$
\int_{-\infty}^{\infty} \delta (t-a)\,f(t)\,dt = f(a),
$$
with $f(t)$ any function. In your case those integrals would thus yield the values one and zero respectively.
For your second question I will only considering linear time invariant systems. In that case the impulse response of such a system can only contains a Dirac delta function if the transfer function of that system has a numerator of the same order as the denominator. Namely, any transfer function of the form
$$
G(s) = \frac{b_n\,s^n + b_{n-1}\,s^{n-1} + \cdots + b_1\,s + b_0}{s^n + a_{n-1}\,s^{n-1} + \cdots + a_1\,s + a_0},
$$
with $b_n \neq 0$ can also be written as
$$
G(s) = b_n + \frac{b'_{n-1}\,s^{n-1} + \cdots + b'_1\,s + b'_0}{s^n + a_{n-1}\,s^{n-1} + \cdots + a_1\,s + a_0},
$$
with $b'_k = b_k - b_n\,a_k$. The inverse Laplace transform of the constant $b_n$ would contribute a Dirac delta term. For the remaining part of the transfer function one could use partial fraction expansion to show that it can't contribute a Dirac delta term.
If a physical system would have a numerator of the same order as the denominator then it would require that the output of the system is directly affected by the input. An example of such physical system might be some electrical motor where you input a voltage and measure the angular position with some voltage leakage from the input signal to the output. However, most physical systems have numerator of a lower order as the denominator. It is more likely that you might encounter equal order numerators and denominators in digital filters (though, those would be z-domain and not s-domain, but roughly the same argument holds) such as notch filters. Those filters are however often used in series with physical systems, so their combined transfer function would also have a lower order numerator.


Image not divisable by 8? (JPEG Compression)

Then it's size is extended by padding zeros or repeating at the end, such that the new size is divisible by 8. Note that reverse repetition at the end reduces the DCT compression size on that block compared to a zero pad approach.

For example if the image was 29 x 27 , then you should add 3 pixels at the right end and 5 lines at the bottom so that it becomes 32 x 32.

Note that you should include the pixel size of the image in the header part so that during decoding process you create the true sized image as 29 x 32


Update Sub Set of the State Vector in Kalman Filter

It will depend upon the measurement model. If your sensor provides data for all the states i.e. for $x_1,x_2,x_3,x_4,x_5$, then you can update only $x_1$ and $x_2$ by making the noise covariance corresponding to $x_3,x_4,x_5$ very large. Actually, there will still be some corrections in these three states but they will be small. The higher the covariance corresponding to these three states, the lower is the correction in them.
Alternatively, you can make the last three rows and last three columns of process covariance very small.


Extract same in phase signal from 3 different sources with shared delayed noise?

Depending on how good your approximation of x(t) is, you might want to try this:

Subtract your approx. of x(t) from all $S$
Cross correlate the signals to get $\phi_2$ and $\phi_3$ (choose $\phi_1$ = 0)
You now should have n(t),
Subtract $n(t)$ from orig $S_1$.
You now have x(t)

edit: Please note, that this will only work if your approximation of x(t) is quite good, and the noise you have really is equal in all mics (and not introduced due to electronics etc.)


How to remove weak edges in a binary photo?

When you are dealing with a binary image (black and white with 1 bit per pixel intensity) you loose all the continuous gradient information, and therefore cannot make use of edge detection (or removal) algortihms that rely on it.

Instead you can return to geometric features such as line length, edge thickness (indeed line thikness), surface area etc to classify features found in the image.

A weak edge should have a corresponding definition in terms of feature geometries so that it can be detected and removed (replaced) if necessary.


I want to know how in practice we can take a kHz/MHz signal to a GHz band?

You guessed right: a mixer.

  Should I build an oscillator with a frequency of 3499700 kHz (or 3480 MHz) to use it as a mixer?

Assuming these frequencies are the difference between your source and your target frequencies, yes.

But you don't use an oscillator as a mixer, you use its output to feed it to a mixer.

  Is it possible in practice to have a precise 3499700 kHz oscillator?

Define "precise"!

Precise is a relative term. If you build a house, your standard of length needs to be accurate to say 0.2 % (so that if your house is 10m long, its length error will be no more than 2cm). If you build an aircraft engine, you might want to be better.

With oscillators it's usually way, way more precise than that, even on cheap standards.

Technologically, RF devices typically use synthesizers to generate local oscillator (LO) frequencies from a reference oscillator, which might be running, for example, at 10 MHz. A synthesizer typically contains all the parts necessary to derive that frequency – a control loop for a voltage-controlled oscillator, that oscillator, some sort of frequency multiplier, filtering to actually get the multiple you want, control and compensation circuitry…

The precision and stability (two different things!) of the resulting oscillator usually mainly depends on the precision and stability of the 10 MHz reference, and to a smaller extent (at least in the frequency regions you mention) on the noise and precision of the control loops.

  Is there any IC useful or should it be by PLL?

This question makes no sense. All the ICs that might make sense here contain at least one form or another of a PLL.

  How?

Good!

  Is there any other and better solution?

To what?

You seem to want a mixer + filters + a synthesizer. Mixers and filters can be bought as modules for moderate prices (e.g. from minicircuits), synthesizers often either take the shape of signal generators in laboratory settings, or as synthesizer ICs + support circuitry in actual devices.


Changing the number of cycles of an sinusoidal input signal

So yes, increasing the duration of the transmitted sinusoidal pulse (by increasing the number of cycles transmitted) will provide improved detection at the same SNR.

This is equivalent to a brute force bit duplication used to decrease the bit error rate (hence improve the reliable communication) at a given SNR. Note, however, that you will be transmitting at a slower rate than a shorter window.

So in effect, this method of longer window will compromise transmission rate vs improved detection at longer distance (lower SNR).


How can I get maximum complex value from an array of complex values(524289) in python?

"Maximum" is ambiguous when it comes to complex values.  The complex value itself doesn't have a min or max.  Which is greater, 1+0j or 0+1j?  What about 1+0j and 0+.5j?  The answer to these questions determines what exactly what you want to do.

You can get the maximum real part (max(1+0j, 0+10j) == 1), maximum imaginary part (max(10+0j, 0+1j) == 1), maximum absolute value (max(10+0j, 0+1j) == 10, max(1+0j, 0+10j) == 10), complex value with the maximum real part (max(1+0j, 0+10j) == 1+0j), complex value with the maximum imaginary part (max(10+0j, 0+1j) == 0+1j), or complex value with the maximum absolute value (max(10+0j, 0+1j) == 10+0j, max(1+0j, 0+10j) == 0+10j).  All are possible with numpy arrays.

The default in numpy (arr.max() or np.max(arr) when arr is complex) is the complex value with the maximum real part.

import numpy as np

arr = np.random.random(1000)+np.random.random(1000)*1j  # generate example data

maxreal = arr.real.max()  # maximum real part
maximag = arr.imag.max()  # maximum imaginary part
maxabs = np.abs(arr).max()  # maximum absolute value

maxcompreal = arr[arr.real.argmax()]  # complex value with maximum real part
maxcomp = arr.max()  # complex value with maximum real part, same as above
maxcompimag = arr[arr.imag.argmax()]  # complex value with maximum imaginary part
maxcompabs = arr[np.abs(arr).argmax()]  # complex value with maximum absolute value


Are all isotropic bandpass filter kernels, wavelets?

Certainly not. ,,Isotropic'' means that the impulse response is independent of direction, that is, that the filter is rotationally symmetric. Wavelets are orthogonal basis functions (filters, maybe) with certain properties. Not every orthogonal basis function will be isotropic, and not every isotropic filter is orthogonal. It is even hard to see how a isotropic filter may be a wavelet.

Look up ,,wavelet conditions'' to learn more about how this orthogonality has to be understood to build an invertible transform. Banded Matrices with Banded Inverses might help the algebraically inclined.


Vector of random numbers with fixed sampling rate

Sampling rate doesn't change anything. Simply call rand function:

x = rand(1, N);


Why samples per symbols $\geq 2$ in the GNURadio Constellation Modulator block?

The GNURadio Constellation modulator enforces a minimum of 2 Samples Per Symbol mostly due to simplicity and some practical reasons. Theoretically, if you're just transmitting PSK without any pulse shaping, you could do 1 sample per symbol (note these are complex symbols) . But typically we like filtered PSK, and the block does this with the RRC filter. However - a key point is that when filtering the signal, you have bandwidth expansion. So if you use a 35% excess bandwidth, your sample rate must be at least 1.35x larger , or another way to say it is that you need at least 1.35 samples per symbol to satisfy nyquist. Similarly for 50% bandwidth you'd need 1.5 samples/symbol. The easiest way to deal with this is to just upsample by 2 (i.e. oversample) and filter with an RRC designed for 2 samples per symbol and then you can support the full range of excess BW values that the RRC can do.
There is one other key idea here though - if you were to try doing a rational upsampling by 1.35x for a 35% excess BW RRC, your resulting signal would be fully occupying your bandwidth (i.e. you have no oversampling). This becomes a problem in multiple areas - for one, when you upsample by 2 later on, the image that must be filtered will show up right next to the one you want to keep, making your digital filtering requirements harder. Additionally, if you instead plan to send this out via a D/A, the image in your 2nd nyquist zone will be sitting right next to your desired image in zone 1, and again filtering with your analog reconstruction filter will be very difficult (e.g. a complex filter with steep transition band).
Also take a look at this very similar question asked by myself where Dan Boschen gave a great answer: Sample rates, Samples per Symbol, and Digital Pulse Shaping


Python equivalent code for Matlab bandpass function

I know it is a bit late, but i had the SAME issue as you today. But from what you showed, i may got the solution similar to Matlab. I don't know if it is the right answer but you could tell me in a future!
Just to complement what you have already said, in the matlab documentation for the function, we have: stopband attenuation of 60 dB.
I did some changes from my variables, to fit your variables.

from scipy.signal import firwin, remez, kaiser_atten, kaiser_beta, kaiserord
Fs = 128.0
lowcut = 1 #In Hertz
highcut = 5 #In Hertz
stopbbanAtt = 60  #stopband attenuation of 60 dB.
width = .5 #This sets the cutoff width in Hertz
nyq = 0.5*Fs
ntaps, gb= kaiserord(stopbbanAtt, width/nyq)
atten = kaiser_atten(ntaps, width/nyq)
beta = kaiser_beta(atten)
a = 1.0
taps = firwin(ntaps, [lowcut, highcut], nyq=nyq, pass_zero=False,window=('kaiser', beta), scale=False)
filtered_signal = signal.filtfilt(taps, a, signal)

From that, i applied FFT in both: Matlab & Python and in spectral analisys i've got the same cut frequencies with my own 1-D data from EEG .


The order is classically the maximum delay (in samples) that is  used to produce each output (Filter order by JOS). For FIR filters, the length is "one plus order":

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

would be 3-tap and order 2.
If the filter is IIR, it would have "more coefficients" than "one plus order". I am not sure one can say "taps" for IIR filter. This recursive  equation would be order 2:

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

A $32^{\mathrm{nd}}$-order with infinite impulse response (IIR) is likely to have stability issues. Even if it does not, the technique that goes "with each tap randomly varied by up to 5% to simulate physical differences" would be likely to endangered the stability randomly.


Determining if the system is linear and time invariant

Hint : Use following properties of Discrete-Time Fourier Transform :

If $DTFT\{x[n]\} = X(e^{j\omega})$, then $DTFT\{x[n - n_0]\} = e^{-j\omega n_0}X(e^{j\omega})$
$IDTFT\{j\frac{\partial{}}{\partial\omega}X(e^{j\omega})\} = n x[n]$


Why does the high frequency component of an image contain more noise?

There are different types of noise sources in imaging process.

Photoelectronic

photon noise: Poisson distribution, whole spectra;

thermal noise: Gaussian distribution, white noise, whole spectra;

Impulse

salt and pepper noise: mainly high frequency noise; (@user7358 made a very good point on the reason).

line drop: part or all of a line in the image lost, can be viewed as a deterministic high frequency noise;

Structured

periodic: low frequency noise;

aperiodic: JPEG noise, low frequency too;

detector striping: calibration differences among individual scanning detectors, mainly low frequency;

If you only observe the evident noise when you zoom in the high frequency band with wavelet, I assume it is mainly the shot noise (salt and pepper). When you are doing some simulations with noise added, try (1) salt and pepper,  (2) Gaussian, and (3) Gaussian with a high-pass band applied, and observe the difference after wavelet decomposition.


In Convolutional Neural Nets, what do convolutions look like?

The kernels used by a ConvNet are nothing but neural weights. You can think of them as a multilayer perceptron with some connections cut off and some weights restricted to be equal (weight sharing).

With this in mind, we must take into account that the kernels (or filters, in this context) are learned, so they depend exclusively on the type of inputs and the desired outputs. Each learning process is unique for a set of inputs/outputs, and so there is no way to know a priori what the filters will turn out to be. They are, as you state in your question, application-dependant.

They can be interpreted, though. Take for example a simple case, where we are trying to classify images that consist of horizontal lines only, and others of vertical lines only. If you do this, then the kernels will indeed look like that: vertical and horizontal lines. That's because the filters have learned common patterns to look for on the inputs, and they have acquired that form in order to be able to match those patterns and maximize their correlation to them.

In Figure 3 at Krizhevsky et al., there are some kernels for you to see. Kernels looking for horizontal lines will consist of horizontal lines. Kernels looking for vertical lines will consist of vertical lines. And so on.

To see this more clearly, you can train a CNN and take one of the kernels learned and filter any image with it to see what happens.

I think that a paragraph from this website can be of relevance:

  So what can we conclude from these feature maps? It's clear there is spatial structure here beyond what we'd expect at random: many of the features have clear sub-regions of light and dark. That shows our network really is learning things related to the spatial structure. However, beyond that, it's difficult to see what these feature detectors are learning. Certainly, we're not learning (say) the Gabor filters which have been used in many traditional approaches to image recognition. In fact, there's now a lot of work on better understanding the features learnt by convolutional networks. If you're interested in following up on that work, I suggest starting with the paper Visualizing and Understanding Convolutional Networks by Matthew Zeiler and Rob Fergus (2013).

In the paper mentioned in the previous quote, there are lots of examples of learned kernels, related to the images that were used to train the network. I think that paper can give you some useful insight.


You reasoning is NOT correct. One period is 10 samples long, not 11. The 11th sample is the first sample of the 2nd period, not the last sample of the 1st period.

So originally you have

sample rate: 10 MHz
sine wave frequency: 1 MHz, which is 1/10th of the sample rate

If you want to get a 1 kHz sine wave you need to set the sample rate to 10 times 1 kHz, i.e. 10 kHz. If you play it back at a sample rate of 909 Hz, your resulting frequency would be 90.9 Hz.


The result looks right, though you've put "1 impulse = ...", when it should be "1 impulse/ms = ..." in the 3rd equation.

Have you access to the equipment to check how fast the fluid actually flows and/or how many milliseconds there are between pulses from the sensor?


Calculate input impedance from S-parameters

S11 is in the frequency domain, and the impedance should also be in the frequency domain. If you want to resample the result, simply use the resample command:

z0           = 50;
A            = importdata("sparameters_lowpass.dat");
N            = 5000;      % number of samples in final result desired

freq         = A(:,1);
realPart     = A(:,2);
imagPart     = A(:,3);

s11_freq     = realPart + imagPart*i; % In frequency domain

z_in         = z0*(1+ (s11_freq))./(1- (s11_freq));

freq_resamp = resample(freq, N)
z_in_resamp = resample(z_in, N)

Note that your formula for S11 would be correct when it is given as a complex signal with a magnitude between 0 and 1 and any phase. Just a word of caution in case the result is in dB magnitude and phase (I don't think that is the case given the extraction of the real and imaginary part shown but that should be confirmed).


Experimental determination of impulse response function

You can determine the impulse response of a system using Gaussian random noise as a test signal, let's name it $x$. Then, you have to measure the system response $y$ to that noise. With a considerable amount of digital samples from both signals then you can compute the crosscorrelation between $x$ and $y$ as follows:

$$
R_{xy}(n\Delta t) = \frac{1}{N-n} \sum_{k=1}^{N-n-1} x[k] \cdot y(k-n\Delta t)
$$

If the white noise has a bandwidth that is much greater than the bandwidth of the system (for instance, 100 times) then the crosscorrelation is equal to the impulse response times the power spectrum density of the noise (which is constant) divided by 2. Of course, there is some uncertainty in the estimation, which decreases with large values for $N$.

Source: Oliver-Cage: Electronic Measurements And Instrumentation


Is it possible to add some kind of noise to a non stationary (having unit root) signals to make them stationary?

One processing approach is to use the Two-Sample Variance (also called the Allan Variance), which takes the differencing over intervals of increasing duration from which it can be directly observed over what interval the original signal (unprocessed) can be assumed to be stationary, and for other durations this processing can convert non-stationary signals to be stationary for further statistical evaluation. This is used widely in the Clock world due to the non-stationary characteristics of oscillator phase noise, but is applicable to other signals as well that have similar non-stationary noise processes (1/f noise, drift, etc).
More information on the Allan Variance (or it’s square root the Allan Deviation, ADEV) is here at Wikipedia:
https://en.wikipedia.org/wiki/Allan_variance
Differencing is a high pass comb filtering response, and works in the case of converting non-stationary clock frequency noise since the non-stationary aspects of the noise are more pronounced as you get closer and closer to the carrier. You can also see this with daily stock market prices where the cost vs time is a random walk process (non-stationary) but the daily difference would be a stationary white Gaussian noise process (assuming no memory in the market which is typically valid for prices from one day to the next- even from one minute to the next which is why traders are concerned with  ns time intervals with regards to trading information / communications)
And I have other posts here on DSP.SE that demonstrate the application and utility of the Allan Variance and Allan Deviation beyond the evaluation of Clocks:
see
How to interpret Allan Deviation plot for gyroscope?
and
What determines the accuracy of the phase result in a DFT bin?


What is meant by a system's "impulse response" and "frequency response?"

The impulse response and frequency response are two attributes that are useful for characterizing linear time-invariant (LTI) systems. They provide two different ways of calculating what an LTI system's output will be for a given input signal. A continuous-time LTI system is usually illustrated like this:

In general, the system $H$ maps its input signal $x(t)$ to a corresponding output signal $y(t)$. There are many types of LTI systems that can have apply very different transformations to the signals that pass through them. But, they all share two key characteristics:

The system is linear, so it obeys the principle of superposition. Stated simply, if you linearly combine two signals and input them to the system, the output is the same linear combination of what the outputs would have been had the signals been passed through individually. That is, if $x_1(t)$ maps to an output of $y_1(t)$ and $x_2(t)$ maps to an output of $y_2(t)$, then for all values of $a_1$ and $a_2$,

$$
H\{a_1 x_1(t) + a_2 x_2(t)\} = a_1 y_1(t) + a_2 y_2(t)
$$

The system is time-invariant, so its characteristics do not change with time. If you add a delay to the input signal, then you simply add the same delay to the output. For an input signal $x(t)$ that maps to an output signal $y(t)$, then for all values of $\tau$,

$$
H\{x(t - \tau)\} = y(t - \tau)
$$

Discrete-time LTI systems have the same properties; the notation is different because of the discrete-versus-continuous difference, but they are a lot alike. These characteristics allow the operation of the system to be straightforwardly characterized using its impulse and frequency responses. They provide two perspectives on the system that can be used in different contexts.

Impulse Response:

The impulse that is referred to in the term impulse response is generally a short-duration time-domain signal. For continuous-time systems, this is the Dirac delta function $\delta(t)$, while for discrete-time systems, the Kronecker delta function $\delta[n]$ is typically used. A system's impulse response (often annotated as $h(t)$ for continuous-time systems or $h[n]$ for discrete-time systems) is defined as the output signal that results when an impulse is applied to the system input.

Why is this useful? It allows us to predict what the system's output will look like in the time domain. Remember the linearity and time-invariance properties mentioned above? If we can decompose the system's input signal into a sum of a bunch of components, then the output is equal to the sum of the system outputs for each of those components. What if we could decompose our input signal into a sum of scaled and time-shifted impulses? Then, the output would be equal to the sum of copies of the impulse response, scaled and time-shifted in the same way.

For discrete-time systems, this is possible, because you can write any signal $x[n]$ as a sum of scaled and time-shifted Kronecker delta functions:

$$
x[n] = \sum_{k=0}^{\infty} x[k] \delta[n - k]
$$

Each term in the sum is an impulse scaled by the value of $x[n]$ at that time instant. What would we get if we passed $x[n]$ through an LTI system to yield $y[n]$? Simple: each scaled and time-delayed impulse that we put in yields a scaled and time-delayed copy of the impulse response at the output. That is:

$$
y[n] = \sum_{k=0}^{\infty} x[k] h[n-k]
$$

where $h[n]$ is the system's impulse response. The above equation is the convolution theorem for discrete-time LTI systems. That is, for any signal $x[n]$ that is input to an LTI system, the system's output $y[n]$ is equal to the discrete convolution of the input signal and the system's impulse response.

For continuous-time systems, the above straightforward decomposition isn't possible in a strict mathematical sense (the Dirac delta has zero width and infinite height), but at an engineering level, it's an approximate, intuitive way of looking at the problem. A similar convolution theorem holds for these systems:

$$
y(t) = \int_{-\infty}^{\infty} x(\tau) h(t - \tau) d\tau
$$

where, again, $h(t)$ is the system's impulse response. There are a number of ways of deriving this relationship (I think you could make a similar argument as above by claiming that Dirac delta functions at all time shifts make up an orthogonal basis for the $L^2$ Hilbert space, noting that you can use the delta function's sifting property to project any function in $L^2$ onto that basis, therefore allowing you to express system outputs in terms of the outputs associated with the basis (i.e. time-shifted impulse responses), but I'm not a licensed mathematician, so I'll leave that aside). One method that relies only upon the aforementioned LTI system properties is shown here.

In summary: For both discrete- and continuous-time systems, the impulse response is useful because it allows us to calculate the output of these systems for any input signal; the output is simply the input signal convolved with the impulse response function.

Frequency response:

An LTI system's frequency response provides a similar function: it allows you to calculate the effect that a system will have on an input signal, except those effects are illustrated in the frequency domain. Recall the definition of the Fourier transform:

$$
X(f) = \int_{-\infty}^{\infty} x(t) e^{-j 2 \pi ft} dt
$$

More importantly for the sake of this illustration, look at its inverse:

$$
x(t) = \int_{-\infty}^{\infty} X(f) e^{j 2 \pi ft} df
$$

In essence, this relation tells us that any time-domain signal $x(t)$ can be broken up into a linear combination of many complex exponential functions at varying frequencies (there is an analogous relationship for discrete-time signals called the discrete-time Fourier transform; I only treat the continuous-time case below for simplicity). For a time-domain signal $x(t)$, the Fourier transform yields a corresponding function $X(f)$ that specifies, for each frequency $f$, the scaling factor to apply to the complex exponential at frequency $f$ in the aforementioned linear combination. These scaling factors are, in general, complex numbers. One way of looking at complex numbers is in amplitude/phase format, that is:

$$
X(f) = A(f) e^{j \phi(f)}
$$

Looking at it this way, then, $x(t)$ can be written as a linear combination of many complex exponential functions, each scaled in amplitude by the function $A(f)$ and shifted in phase by the function $\phi(f)$. This lines up well with the LTI system properties that we discussed previously; if we can decompose our input signal $x(t)$ into a linear combination of a bunch of complex exponential functions, then we can write the output of the system as the same linear combination of the system response to those complex exponential functions.

Here's where it gets better: exponential functions are the eigenfunctions of linear time-invariant systems. The idea is, similar to eigenvectors in linear algebra, if you put an exponential function into an LTI system, you get the same exponential function out, scaled by a (generally complex) value. This has the effect of changing the amplitude and phase of the exponential function that you put in.

This is immensely useful when combined with the Fourier-transform-based decomposition discussed above. As we said before, we can write any signal $x(t)$ as a linear combination of many complex exponential functions at varying frequencies. If we pass $x(t)$ into an LTI system, then (because those exponentials are eigenfunctions of the system), the output contains complex exponentials at the same frequencies, only scaled in amplitude and shifted in phase. These effects on the exponentials' amplitudes and phases, as a function of frequency, is the system's frequency response. That is, for an input signal with Fourier transform $X(f)$ passed into system $H$ to yield an output with a Fourier transform $Y(f)$,

$$
Y(f) = H(f) X(f) = A(f) e^{j \phi(f)} X(f)
$$

In summary: So, if we know a system's frequency response $H(f)$ and the Fourier transform of the signal that we put into it $X(f)$, then it is straightforward to calculate the Fourier transform of the system's output; it is merely the product of the frequency response and the input signal's transform. For each complex exponential frequency that is present in the spectrum $X(f)$, the system has the effect of scaling that exponential in amplitude by $A(f)$ and shifting the exponential in phase by $\phi(f)$ radians.

Bringing them together:

An LTI system's impulse response and frequency response are intimately related. The frequency response is simply the Fourier transform of the system's impulse response (to see why this relation holds, see the answers to this other question). So, for a continuous-time system:

$$
H(f) = \int_{-\infty}^{\infty} h(t) e^{-j 2 \pi ft} dt
$$

So, given either a system's impulse response or its frequency response, you can calculate the other. Either one is sufficient to fully characterize the behavior of the system; the impulse response is useful when operating in the time domain and the frequency response is useful when analyzing behavior in the frequency domain.


What does the amplitude in a frequency response of a signal signify?

Short answer, it depends what exactly you're looking at.

If you have a time series and compute the DFT, this will give you the concentration of energy at each discrete frequency point you have chosen. When you compute the PSD of a signal in MATLAB (using pwelch or cpsd), you'll get units of power/frequency. This describes the amount of energy or power in that specific frequency range (bin). If instead your goal is to relate the amplitudes of sinusoids in a time series directly to the values that you're reading off the peaks in the PSD plot, you need to ensure that your FFT is appropriately scaled. In the MATLAB functions cpsd() or pwelch(), the FFT is scaled such that the $\textit{total energy}$ in the signal is preserved (done by dividing by the FFT by the sampling period, $t_s$). If you want to recover the original amplitudes of a signal (e.g., $x=A_1 sin(\omega_1*pi*t)+A_2 sin(\omega_2*pi*t))$, you would want to ensure that the FFT is scaled by the duration of the signal, $T$.

Basically, to make any accurate inference from the frequency response plot, you need to know exactly how it was computed (units).

Some refs. that I found helpful:
What is the difference between the PSD and the Power Spectrum?


How to undo dynamic Doppler effect by software?

If the only reason for the buffer is so your algorithm has something to process, then the simplest way to implement this is to keep your current algorithm but use a shorter buffer. If the larger buffer was present for other reasons (i.e. data arrives in large chunks), then keep that larger buffer, but feed it into a shorter buffer that is given to your processing code. (It seems any algorithm able to recognize Doppler shift in a given snapshot must also be able to process that same snapshot.) In this way you can keep your current algorithm, and tailor the short buffer length as desired. I can't imagine how an algorithm specifically designed to do this per sample would be of benefit, unless it takes less CPU--which seems unlikely.


Are noncausal filters ever used in practice?

Yes.

The problem with a system that operates in (near) real time is that you can't look into the future. One way you can deal with this is if you only need a finite amount of look-ahead is to put some delays and then delay the output so you're still causal.

However, many filtering problems have non-causality allowed, e.g. filtering a file on a disk (which occurs a lot. Audio or image or video files you download, time series such as finance data, or histories of systems). For example, you can collect a time series and smooth it with the rauch-tung-striebel filter rather than running a fixed lag smoother.


Continuous Wavelet Transform time vector in python

Similar question here Which time-frequency coefficients does the Wavelet transform compute?

The picture you've shown is used for DWT such as pywt.wavedec, not CWT.

CWT is a continuous function; it exists at all points in the time-scale plane.

pywt.cwt produces a 2D array with constant numbers of columns and rows, so it's a regular sampling of the CWT, at whatever scales you specify, for every sample of the input:

shape(y)
Out[5]: (512,)

shape(freqs)
Out[6]: (128,)

shape(coef)
Out[7]: (128, 512)

I think the DWT's coefficients contain just enough information to reproduce the original exactly with an inverse DWT, while the CWT output has a lot of redundant information and is maybe missing some information.


White noise and entropy confusion about concepts

A signal's spectral entropy is not really a measure of information. Rather, it is a statistical measure of dispersion in the signal's power spectrum.

In other words, a random noise signal will tend to have large spectral entropy because its power spectrum is distributed among all frequencies. Speech signals tend to concentrate their power in a small number of bins, and their power spectrums are not flat. This results in lower spectral entropy.

A few comments about entropy:

Signals have information; sources have entropy. Signals never have entropy, and any reference to such concept is an abuse of terminology.
In communications, we prefer to have sources with large entropy, because then the messages that we transmit will convey a large amount of information. If you have a source with low entropy, then you should apply a source encoder to turn it into a large entropy source.
Shannon's entropy deals with uncorrelated sources. Many real-world sources produce correlated signals (such as text or speech) and their entropy is harder to define and calculate.


Why does this constant-Q transform "lose octaves"?

The lower frequency atoms are longer in time (or of a length roughly proportional to the period of the frequency, give or take 2X).  So less of these lower frequency atoms are needed per second to cover that time period.  And more of the shorter atoms (for shorter wavelength frequencies) are required for about the same coverage (or with roughly the same percentage overlap).

And to maintain a nearly constant Q, the length of a bandpass FIR filter ends up roughly proportional to some constant multiple of the wavelength of the center bandpass frequency of that filter.


Estimating number of sinusoid from spectral estimation

The best approach and the details of how to do it depend greatly on the particulars of your situation.  If all your vectors are of the same length, and all your sinusoids have an integer number of cycles per vector (and therefore their frequencies correspond to a DFT bin) you can take the DFT of each vector.  Add the vectors coherently if the phases are the same in each random realization, add the DFT magnitudes if the phases are not the same.  Then set a threshold and find the points above the threshold.

If the frequencies are not on DFT bins but are spaced apart from each other, you need to modify the above approach to take account of the fact that the energy of each sinusoid is spread over multiple DFT bins.  You can still use threshold-based detection but will need to interpolate to find the frequencies.

If the frequencies could be spaced very closely together - which, from your plot, appears to be the case - the DFT-based approach is difficult.  The best approach may be a subspace-based technique like ESPRIT (Estimation of Signal Parameters Via Rotational Invariance Technique).  This involves estimating the covariance matrix and computing its eigendecomposition, then using the eigenvalues to estimate the power of each sinusoid and the eigenvectors to estimate the frequencies.  See, for examples: Richard Roy, Thomas Kailath: "Estimation of Signal Parameters via Rotational Invariance Techniques." IEEE Transactions on Acoustics, Speech, and Signal Processing.  Vol 37, No 7, July 1989.


How to design an all pass type filter with inverted phase response?

Sort of. The inverse of an allpass filter is also an allpass filter and can simply be calculated by time flipping the impulse response. However that makes the filter non-causal. So you need to add bulk delay to make it causal again. As long as you don't need to do real time processing, that's typically not a problem.

EDIT: add a specific example. The key is to do the filtering in the time domain. You need to truncate the IIR response to an FIR response using "enough" samples.

% create an arbitrary allpass filter
b = [-.5 1]; a = [1 -.5];

% calculate the impulse response
n = 128; d0 = zeros(n,1); d0(1)= 1;
h1 = filter(b,a,d0);

% time flip and shift by 1 to keep the first sample in the same spot
h2 = circshift(flipud(h1),1);

% plot frequency responses
fh1 = fft(h1);
fh2 = fft(h2);
phases = unwrap([angle(fh1(1:n/2+1,:)) angle(fh2(1:n/2+1,:))]);
plot(phases);


How to statistically quantify the coverage of features in the data?

To visualize the features, why don't you just do scatter plots? I assume that if you do feature extraction, you know where in the image the features are found. If different features have different "sizes", i.e. some features are point features, whereas others take up e.g. a 5x5 window (and are non-overlapping), you may want to adjust the size of your symbols for the scatter plot to visually improve the relative amount of coverage.

If you want to quantify the coverage, then you're looking at two independent measures

Density: To this end, you can simply divide the number of features by the image size.
Homogeneity: I've seen Ripley's K used to great success when looking at whether data is homogeneously distributed, or whether it is clustered.


Expected value of coupled signals is zero

It means u and v are orthogonal.

If u or v or both have zero mean then they are also uncorrelated.

If they are uncorrelated and jointly gaussian then they are also independent.


Identify the system from input and output

You can use an equation error method and get an estimation of the coefficients by solving an overdetermined system of linear equations. The transfer function of an IIR filter is given by

$$H(z)=\frac{B(z)}{A(z)}=\frac{\sum_{m=0}^{M}b[m]z^{-m}}{1+\sum_{n=1}^{N}a[n]z^{-n}}\tag{1}$$

From (1), for given input and output sequences $x[n]$ and $y[n]$, respectively, the following equation must be satisfied:

$$(a*y)[n]=(x*b)[n]\tag{2}$$

where $*$ denotes convolution. Equation (2) can be written in matrix/vector form:

$$\mathbf{y}=\mathbf{Xb}-\mathbf{Ya}\tag{3}$$

where $\mathbf{y}$ is the vector of the output signal, $\mathbf{b}$ and $\mathbf{a}$ are the unknown coefficient vectors (excluding $a[0]=1$), and $\mathbf{X}$ and $\mathbf{Y}$ are the convolution matrices consisting of the values of the input and output sequences, respectively. Equation (3) can be rewritten as

$$\mathbf{y}=\mathbf{[-Y,X]\cdot[a,b]^T}\tag{4}$$

where $^T$ denotes transposition. The overdetermined system (4) can be solved in a least squares sense using standard methods, such as the one implemented in Matlab or Octave.

Note that in any case you need to estimate/choose the orders $M$ and $N$ of the numerator and denominator polynomials in (1). For most standard filters you have $M=N$.

This little Matlab/Octave script shows how it works:

x=randn(100,1);       % some input signal
a=[1 .3 -.4 .2 .1];   % denominator coefficients
b=[2 3 4 3 2];        % numerator coefficients
y=filter(b,a,x);      % filter output
M=3;                  % guessed numerator order
N=3;                  % guessed denominator order

% build matrices and solve overdetermined linear system
A = flipud(hankel(flipud(y(:)))); A1 = A(:,1); A = A(:,2:N+1);
B = flipud(hankel(flipud(x(:)))); B = B(:,1:M+1);
C = [-A,B];
c = C\A1;
a1 = [1;c(1:N)];
b1 = c(N+1:M+N+1);

% compare outputs of original and estimated filter
y1=filter(b1,a1,x);
n=1:100;
plot(n,y,n,y1)


What does the frequency band mean when it comes to finding aliases?

Let's use a simple example: A sine wave of 400 Hz. The Fourier spectrum has two components: one at 400 Hz and one at -400 Hz (since it's a real valued time signal).

Now let's sample this at 1 kHz. Sampling in one domain is periodic repetition in the other domain. So the sampled signal consists of your two original frequencies, -400 Hz and +400 Hz repeated again and again in 1 kHz intervals. So you get components at +-400 Hz, +-600 Hz, +-1400 Hz, +-1600 Hz, etc.

Any of these mirror frequencies (600, 1400, 1600, 2400, 2600, ...) will give you the exact same samples as the original 400 Hz.


Noise problem when transforming and quantizing difference between images

JPEG was developed for photographic content, not for difference-of-pictures content. Therefore, the fact that it's not overly applicable to the difference between successive frames isn't surprising.

In your specific quantization approach, where you eradicate high frequencies, that's actually even pretty intuitive: The difference between successive frames are probably small, and fine-grained. Of course, when you erase high frequency content, you lose most of the difference information and get a large error. So, I'd call this an indication that your high-frequency deletion works, and little more.

So, wrong approach to video coding; you could use JPEG for I-frames and code the differences to previous frames with some sort of linear predictor. Or, look at existing video codecs, and adapt some of the methods used there - in the case of grass, which presumably is blowing in the wind in a pretty uniform direction with smaller local variations, motion prediction coding might be of special interest to you.


How to solve this question by using convolution?

Assuming you mean that $y[n] = h[n] * x[n]$, then
\begin{equation}
\begin{split}
 y[n]  &\stackrel{(a)}{=} h[n] * x[n]\\
 &\stackrel{(b)}{=} \sum_{m=-\infty}^{\infty}
 h[m]x[n-m]\\
 &\stackrel{(c)}{=} \sum_{m=-\infty}^{\infty}
 u[m]2^{n-m}u[n-m]\\
 &\stackrel{(d)}{=}2^n
 \sum_{m=-\infty}^{\infty}
 u[m]2^{-m}u[n-m]\\
 &\stackrel{(e)}{=}2^n
 \sum_{m=0}^{\infty}
 2^{-m}u[n-m]\\
 &\stackrel{(f)}{=}2^n
 \sum_{m=0}^{n}
 2^{-m}\\
 &\stackrel{(g)}{=}2^n
 \frac{1 - (\frac{1}{2})^{n+1}}{1 - \frac{1}{2}}\\
 &\stackrel{(h)}{=}2^{n+1} - 1
\end{split}
\end{equation}

(a) is the convolution symbol.

(b) is using the definition of discrete convolution.

(c) replacing the quantities.

(d) $2^n$ is independent of the summation index.

(e) $u[m]=1$ if $m\geq 0$ and $0$ otherwise.

(f) if $n-m<0$, or if $n<m$ then the summation is zero. So we have (for non-zero quantities), $m = 0 \ldots n$.

(g) Using geometric series

(h) rearrangements.


Effect of Windowing Function on Autocorrelation Function

When you are using the FFT to get a measurement of the frequency content of a signal it is appropriate to window the data first to avoid distortions due to a rectangular cutoff.

That is not what you are doing with the FFT, though.  You are using the FFT and inverse FFT to efficiently calculate a correlation (an autocorrelation is just a special case of correlation).  Done properly with padding and such, the FFT will get you the exact same results that you would get by calculating the correlation in the time domain.  If you window, though, you will not get the same results.  In short, you should not window when you are using the FFT to calculate correlations or convolutions.


The problem occurred due to the pixel values of temp had the type uchar.

The solution to my problem was to output the pixel values as uchar, manipulate them as double, and store them as uchar.

double den=0,num=0;
      for(int i = -(kernel/2); i <= (kernel/2) ; i++)
      {
        for(int j = -(kernel/2) ; j <= (kernel/2); j++)
        {
            den += pow(temp.at<uchar>(row+i,col+j),P);
            num += pow(temp.at<uchar>(row+i,col+j),P+1);
        }
      }

    dst.at<uchar>(row,col) = num/den;


Using a low pass filter to interpolate signal

Another view which results in superior interpolation filter design is revealed by reviewing the spectrum of the signal with zeros inserted (resulting in an increase of the sampling rate by $I$ when you insert $I-1$ zeroes) and you will see that the original spectrum is replicated at integer multiples of the original sampling rate. The ideal low pass filter is the one that can pass the original spectrum with no distortion and reject all the new alias copies perfectly. (Ideal is not achievable but this defines the filter design targets to approach that).  The ideal filter will interpolate the non zero values to grow the zeros to the value for the final interpolated waveform.

What is occurring, to explain the OP's question conceptually, is that the signal requires high frequency content when going from zeros to a large value in just one sample (and inserting zeros causes this to occur as evidenced by the additional high frequency content in the spectrum). This makes complete sense since the frequency content would be related to a change in magnitude versus a change in time (a big change in a short time requires high frequencies). A filter that removes the high frequencies results in the inability for the signal to change rapidly from one sample to the next (due to the memory of the previous samples which is what a filter does). A very easy way to see this occurring is with a simple moving average filter over $M$ samples when you interpolate by $M$ by inserting $M-1$ zeros; such a filter will hold the last sample as a zero-order hold, growing each zero to the last non-zero sample. This isn't a good way to do interpolation since only one non-zero sample is in the memory of the filter, but very easy to visualize. The "magic" happens with better designed filters as described in this answer which then consider many more samples to perform what is essentially interpolation with higher order polynomials.

With regards to Sinc reconstruction that is commonly used to explain "perfect" reconstruction; it is worth mentioning that Sinc reconstruction is only perfect when you can use a Sinc, but unfortunately the Sinc function extends to infinity which is impossible to do in practical realization (it is equivalent to filtering with a brickwall filter). Therefore the Sinc response is truncated in practice, which on its own results in a very poor filter, which is then typically windowed to improve performance. Ultimately, its use when truncated and windowed results in inferior reconstruction filters. Given the approach of paying attention to where the aliases actually reside in frequency, you can realize optimized filters (superior to simply truncating and windowing a Sinc function) using the multi-band filter algorithms such as least squares using firls() function that is available in MATLAB, Python and Octave.  I demonstrate this in the plots below from an X4 interpolation example with $f_s = 10$ KHz. The upper plot compares options for interpolation filters with the same number of taps, showing the poor performance of a truncated Sinc in green, and then the improved windowed Sinc filter in green along with the least squares multiband filter in red. All have the same number of taps, but the multiband approach results in 10 dB better rejection of the higher frequency images (and not visible from the scale of the plot, but lower passband ripple distortion as well).

The filter magnitude for the multiband filter was shifted down to overlay with the spectrum after zero-insert showing how the multiband filter maximizes rejection where needed thus resulting in the least amount of distortion in the resulting interpolated signal.

So the spectrum above with it's replication at the higher frequencies (images) has the zeros inserted in the time domain. When those higher frequencies are removed through filtering, the spectrum represents the original signal, at the higher sampling rate, as if it was sampled directly.

A truncated and properly windowed Sinc filter will not achieve the rejection or passband performance of this filter with the same number of taps since it would provide rejection at all locations in between the aliases, while the multi-band filter maximizes the rejection where it is needed. The resulting distortion is predictable and can be traded with the desired filter complexity. This is a standard interpolation design approach for high performance and efficient interpolation (especially when converted to polyphase structures!).


Are the components of a Bandpass signal namely "inphase" and "quadrature" real functions in time domain?

$X(t) = Re\{XX(t)e^{j2\pi ft}\}$. So $XX(t)$ can be any generic complex baseband signal with real I and Q components - $XX(t) = XX_i(t)+j XX_q(t)$ and $e^{j2\pi ft} = \cos(2\pi ft) + j\sin(2\pi ft)$

After the complex multiplication, and taking real part, you get $X(t) = XX_i(t)\cos(2\pi ft) - XX_q(t)\sin(2\pi ft)$.

So $XX_i(t)$ and $XX_q(t)$ are indeed real functions and their spectra is even-symmetric.

If $\hat{X}(t)$ is the Hilbert transform of $X(t)$, then $X^+(t) = X(t)+j\hat{X}(t)$. $XX_i(t)=X(t)\cos(2\pi ft)+\hat{X}(t)\sin(2\pi ft); XX_q(t)=\hat{X}(t)\cos(2\pi ft)-X(t)\sin(2\pi ft)$


What effect does a delay in the time domain have in the frequency domain?

The discrete Fourier transform (DFT), commonly implemented by the fast Fourier transform (FFT), maps a finite-length sequence of discrete time-domain samples into an equal-length sequence of frequency-domain samples. The samples in the frequency domain are in general complex numbers; they represent coefficients that can be used in a weighted sum of complex exponential functions in the time domain to reconstruct the original time-domain signal.

These complex numbers represent an amplitude and phase that is associated with each exponential function. Thus, each number in the FFT output sequence can be interpreted as:

$$
X[k] = \sum_{n=0}^{N-1} x[n] e^{\frac{-j 2 \pi n k}{N}} = A_k e^{j \phi_k}
$$

You can interpret this as follows: if you want to reconstruct x[n], the signal that you started with, you can take a bunch of complex exponential functions $e^{\frac{j 2 \pi n k}{N}}, k = 0, 1, \ldots, N-1$, weight each one by $X[k] = A_k e^{j \phi_k}$, and sum them. The result is exactly equal (within numerical precision) to $x[n]$. This is just a word-based definition of the inverse DFT.

So, speaking to your question, the various flavors of the Fourier transform have the property that a delay in the time domain maps to a phase shift in the frequency domain. For the DFT, this property is:

$$
x[n] \leftrightarrow X[k]
$$
$$
x[n-D] \leftrightarrow e^{\frac{-j2 \pi k D}{N}}X[k]
$$

That is, if you delay your input signal by $D$ samples, then each complex value in the FFT of the signal is multiplied by the constant $e^{\frac{-j2 \pi k D}{N}}$. It's common for people to not realize that the outputs of the DFT/FFT are complex values, because they are often visualized as magnitudes only (or sometimes as magnitude and phase).

Edit: I want to point out that there are some subtleties to this rule for the DFT due to its finiteness in time coverage. Specifically, the shift in your signal must be circular for the relation to hold; that is, when you delay $x[n]$ by $D$ samples, you need to wrap the last $D$ samples that were at the end of $x[n]$ to the front of the delayed signal. This wouldn't really match what you would see in a real situation where the signal just doesn't start until after the beginning of the DFT aperture (and is preceded by zeros, for example). You can always get around this by zero-padding the original signal $x[n]$ so that when you delay by $D$ samples, you just wrap around zeros to the front anyway. This relationship only applies to the DFT since it is finite in time; it's does not apply to the classic Fourier transform or discrete-time Fourier transform.


Simple way to replicate bloom (glow) effect in image using Matlab?

Work in a linear gamma (intensity-linear) color space. Create a Gaussian blurred version of the image. I don't have MATLAB to try it, but it can be done using imgaussfilt. Then add some of the blurred image to the original image, or mix the two images to taste.

Figure 1. Bloom simulated in Gimp by adding to the image its Gaussian blurred copy at 31.5 % opacity.