What are the advantages and disadvantages separability of Gaussian 2D filter?

Assuming the simplest case with a square image $x[n,m]$ of size $N \times N$ and a square filter kernel $h[n,m]$ of size $M \times M$, the raw 2D convolution to produce the, cropped, output image $y[n,m]$ of size $N \times N$ requires about $N^2 \times M^2$ MACs.

The raw 2D convolution between $x[n,m]$ and the filter $h[n,m]$ is implemented by using two for loops to range through each output pixel and two additional for loops to perform the 2D convolution at that pixel location. Hence a total of 4 nested for-loops are required; a complexity of $\mathcal{O}(N^2M^2)$

When the filter kernel $h[n,m]$ is separable, such as $h[n,m] = f[n]g[m]$, then in this case the convolution between the image $x[n,m]$ and the filter $h[n,m]$ can be performed without a raw 2D convolution sum, by the following approach:

First, perform a 1D convolution between columns of $x[n,:]$ and the 1D filter $f[n]$, which requires about $N \times M$ MACs to complete. This operation should be performed for each column of $x[n,m]$ by proceeding along its horizontal, $N$ many, columns. Hence a total of $N \times M \times N$ MACs will be required to complete the first step to produce the intermediate image $w[n,m]$
Then, apply the similar algorithm, with the filter $g[m]$ and rows of the intermediate image $w[:,m]$, which will require similar number of MACs as $N \times M \times N$.

Hence as a total you will need about $2N^2M$ MACs for the separable implementation of the 2D convolution. The actual number depends on the cropping type you apply. Thus a complexity of about $\mathcal{O}(N^2M)$ is attained.

It can be seen that, to implement the separable convolution algorithm, you will only need 3 nested for loops, which is where the gain comes from.

The main advantage of separable filtering is quite clear; much reduced computational cost. In fact even the 2D-FFT algorithm makes use of it as the 2D-DFT kernel is separable.

Another advantage is that since there are less number of MACs involved to produce the same output compared to the raw 2D convolution, it's less prone to numerical issues and likely to produce more accurate results.

A disadvantage of separation is that it requires an extra ram memory to store the intermediate image, which can be a concern in certain applications.


How to "scale" the FFT when using it to calculate discrete convolution?

this formula $FFT^{-1} \left( FFT \left( F \right) \odot FFT \left( G \right) \right) \stackrel{!}{=} F \ast G$ is not always true, it's true when you scale at the inverse operation as you mentioned in your first case, in second case you have to multiply the right side by $\sqrt{N}$ to make it right.

Generally when you apply the Fourier transform over the convolution you see any scaling would be used for the first signal and no scaling remain for the second signal so if you multiply and divide by the scaling and use the scaling in the denominator for standard Fourier transform, the scaling remain in the nominator.


You have aliasing for sampling signals of any frequency: a signal of 5Hz sampled at 20Hz is indistinguishable from a signal of 15Hz of opposite phase.

It's common in audio usage to work with "lowpass signals" only and consider "aliasing" for frequencies more distant from DC, but in reality sampling of anything is ambiguous by multiples of the sample frequency regardless of the original signal.  You can perfectly well sample signals of 15Hz with a sampling frequency of 20Hz as long as you make sure that the signal only contains frequencies between 10Hz and 20Hz (for example).

For any detected signal component, increasing/decreasing its frequency by multiples of the sampling frequency will not cause any change at the sampled time points.  Hence you want to prefilter the signal such that adding multiples of the sampling frequency to any signal component leads to something that cannot be present in the signal: that way the sampled values represent the signal uniquely.


can't we use other function such as: $R=1-(\lambda_1+\lambda_2)^2$

This formula isn't as useful, because it doesn't differentiate between two equally large eigenvalues (a corner) and one twice as large eigenvalue plus one very small one (a straight line). So for example $\lambda_1 = 2, \lambda_2 = 0$ will give the same result as $\lambda_1 = 1, \lambda_2 = 1$, although the former is probably not a corner.

But generally speaking: Yes, the Harris formula isn't the only one. You can for example use the harmonic mean between the eigenvalues, i.e. $2/(1/\lambda_1 + 1/\lambda_2) = 2 det(M)/trace(M)$. Or you can use the smaller eigenvalue directly.

The Harris corner detector was developed at a time when floating point division and square root were still expensive operations. I'd guess that was one reason why a detector that only uses multiplications was preferred.

These days, on a PC, I'm quite sure you can divide floating point numbers faster than you can get them to the CPU or GPU, so it makes sense to consider other corner detectors. On other hardware, like mobile devices or smart cameras the lower cost can still be an important advantage of the Harris corner detector.


The difference is in the scaling of the power spectrum. I suppose you will find that the difference in scaling always equals $L/F_s$, which for the numbers in your question is indeed $8000/500=16$.

Take a look at Power Spectral Density Estimates Using FFT for the correct scaling. If you normalize the FFT result by the FFT length, you need to scale the squared magnitude of the normalized FFT by $L/F_s$. Furthermore, you need a factor of $2$ if you throw away the negative frequencies (I see you did that). Note, however, that you shouldn't apply a factor of $2$ for the bins at DC and at Nyquist, because they are not mirrored at negative frequencies.

Here is a simple Matlab code from the above quoted mathworks page for computing a periodogram-based one-sided power spectrum estimate using the FFT (my comments):

Fs = 1000;                          % sampling frequency (Hz)
N = length(x);                      % even! (=> bin N/2+1 is Nyquist)
xdft = fft(x);
xdft = xdft(1:N/2+1);               % DC to Nyquist (one-sided)
psdx = (1/(Fs*N)) * abs(xdft).^2;   % periodogram scaling
psdx(2:end-1) = 2*psdx(2:end-1);    % scaling for one-sided PSD
                                    % note: no factor 2 at DC and Nyquist


Recursive Implementation of the Gaussian Filter

The answer was simple, the article uses the coefficients value on one hand where the MATLAB implementation on the other.
Namely, a minus sign should be added.

Here's the correct code:

qFactor = 5;

b0Coeff = 1.57825 + (2.44413 * qFactor) + (1.4281 * qFactor * qFactor) + (0.422205 * qFactor * qFactor * qFactor);
b1Coeff = (2.44413 * qFactor) + (2.85619 * qFactor * qFactor) + (1.26661 * qFactor * qFactor * qFactor);
b2Coeff = (-1.4281 * qFactor * qFactor) + (-1.26661 * qFactor * qFactor * qFactor);
b3Coeff = 0.422205 * qFactor * qFactor * qFactor;

normalizationCoeff = 1 - ((b1Coeff + b2Coeff + b3Coeff) / b0Coeff);

vDenCoeff = [b0Coeff, -b1Coeff, -b2Coeff, -b3Coeff] / b0Coeff;

vXSignal = zeros(61, 1);
vXSignal(31) = 10;

vYSignal = filter(normalizationCoeff, vDenCoeff, vXSignal);
vYSignal = filter(normalizationCoeff, vDenCoeff, vYSignal(end:-1:1));

figure();
plot(vYSignal);


Transition table / matrix for convolutional encoder K=3 understandings

The OP's question reveals many mis-understandings of the issues.
The convolutional code in question is a nonsystematic $[2,1]$ convolutional code with constraint length $K = 3$. The specification is incomplete but the complete specification (we also need to know what the generator polynomials are) can be deduced from the state transition diagram that the OP has provided. This particular convolutional code is used as the first example in almost all texts on coding, and so the code is sometimes referred to as the ad nauseam code. Let's pick apart what things mean.  There is 1 input stream of bits entering the encoder during each clock cycle and 2 streams of bits leaving the encoder (in parallel, on two separate wires) during the clock cycle. The encoder itself contains one shift register (not $3$ as the OP thinks) of length K-1 = 2 meaning that it is capable of storing two bits. The encoder also contains XOR gates for computing the output bits from the encoder input bit and the shift register contents. The bits stored in the shift register are the immediate past input bit and the input bit immediately previous to the immediate past input bit, that is, the input bit from two clock cycles ago. Thus, this single two-bit shift register can be in one of $2^{K-1} = 4$ states rather than three shift registers somehow managing to be one of only $4$ states which has puzzled the OP (and would puzzle anyone else too). At the end of the clock cycle, the current bit enters into the shift register, becoming the "immediately past input bit" for the next clock cycle while the bit that served as the  "immediately past input bit" for the current clock cycle moves over one place to become the "input bit two cycles ago" for the next clock cycle.  At the start, the shift register is initialized with two $0$s in it. That is, the encoder is in state $00$ at the beginning of the encoding process.
During any clock cycle, K = 3 bits  viz. the current input bit, the immediately past input bit and the input bit from two clock cycles ago, are available for computing the 2 output bits,.  The state transition diagram of the encoder is a directed graph whose vertices labeled with the state that the encoder is in, and the label of the directed edge from one state to another state denotes "current input/output bits". For example, if the current state is $00$ and the input is a $0$, the encoder outputs $00$ and remains in state $00$. The edge showing this transition is thus labeled $0/00$ and the edge leaves and enters state $00$. On the other hand, if the input is a $1$, the output is $11$ and the encoder moves from state $00$ to state $01$. This the directed edge from $00$ to $01$ is labelled $1/11$. We also deduce from this that whomsoever drew the logic map for the encoder had the input entering the encoder from the right and the shift register shifting leftwards. If the input bits had been entering from the left and the shift register shifting rightwards, the next state would have been labelled $10$; you pays your money and you makes your choice.
Note that there is no edge from state $00$ to state $11$ or vice versa; there is no way that this particular encoder shift register can make a transition from state $00$ to state $11$ or vice versa in one clock cycle. Note also the symmetries in the state transition diagram. The two edges leaving each state have complementary input labels (obviously, if you think about it a bit, pun intended) and complementary output labels too. The two edges entering each state also have complementary output labels (but the same input label; can you see why?).
So, how does the encoder actually work? Well, if the input data stream is
$$ \cdots d[n-4]\quad d[n-3] \quad d[n-2] \quad d[n-1] \quad d[n] \quad d[n+1] \quad \cdots$$
then at the $n$-th clock cycle, the encoder shift register contains $d[n-2]$ and $d[n-1]$ as indicated below via the box, and $d[n]$ is the current input to the encoder shift register.
$$ \cdots d[n-4]\quad d[n-3] \quad\bbox[yellow,5px,border:2px solid red] {d[n-2] \quad d[n-1]} \quad d[n] \quad d[n+1] \quad \cdots$$
During this clock cycle, the encoder produces two code bits
\begin{align}
c_i[n] &= d[n]\oplus d[n-1] \oplus d[n-2],\\
c_q[n] &= d[n]  \oplus d[n-2].
\end{align}
The shift register then shifts, and during the $(n+1)$-th clock cycle the register contents are as shown below and the input data bits is $d[n+1]$.
$$ \cdots d[n-4]\quad d[n-3] d[n-2] \quad \bbox[yellow,5px,border:2px solid red] {d[n-1] \quad d[n]} \quad d[n+1] \quad \cdots$$
The generator matrix for this code is
$${\bf G} = \big[ 1 + x + x^2 \quad 1 + x^2\big]$$
meaning that if we regard the input data stream as a series
$$d(x) = d_0 + d_1x + d_2x^2 + \cdots,$$ then the encoder output streams can  also be regarded as series
\begin{align} c_I(x) &= (1 + x + x^2)d(x)\\
&= d_0 + (d_1 + d_0)x + (d_2 + d_1 + d_0)x^2 + \cdots\\
& \qquad \quad \cdots + (d_n + d_{n-1} + d_{n-2})x^n + \cdots\\
c_Q(x) &= (1 + x^2)d(x)\\
&= d_0 + (d_1)x + (d_2 + d_0)x^2 + \cdots\\
& \qquad \quad \cdots + (d_n + d_{n-2})x^n + \cdots
\end{align}
or, in short,
$$d(x){\bf G} = \big[c_I(x) \quad c_Q(x)\big] = \big[(1+x+x^2)d(x) \quad (1+x^2)d(x)\big].$$
Dedicated DSPissers will replace $x$ by $z^{-1}$ and say "Hey, you are telling me that $c_I(z) = (1+z^{-1} + z^{-2})d(z)$ and $c_Q(z) = (1+z^{-2})d(z)$, and I know that multiplication in the $z$ domain corresponds to convolution in the time domain, and so why don't you name these codes as convolutional codes and save a lot of explanations on dsp.SE?" and that's exactly what it is: the encoder output bits are obtained by convolving the encoder input sequence with the generator polynomial sequences.

I've seen on the internet something weird that the transition table for convolutional encoder $K=3$
is: transitionMatrix = [0 1 1 55; 0 0 1 1; 55 0 55 1; 55 0 55 1]; it was written that value $55$ is value of STATE HISTORY .
I didn't understand this matrix and how this matrix describe the state machine that I attached the photo above for convolutional encoder $K=3$

I too have no idea what that gobbledygook means.
In a comment, the OP asks about implementing a transition matrix. There may be a language barrier here in communicating with each other. I don't know what exactly implementing a transition matrix means or why one would want to go about implementing one. I know about people implementing the encoder from the state transition diagram, whether in hardware with flip-flops and gates, or via FPGAs, or in MATLAB or Python or C++ etc. When the generator matrix is not known, it takes a little analysis of the state transition diagram to figure out what the generator polynomials are, but once that is done, implementing the encoder is quite straightforward. Another way of looking at it is that the encoder implementation implicitly implements the transition matrix in that the transitions that the encoder makes follows the rules laid out by the transition matrix.


Octave is a multi-platform open source math and matrix toolkit. It has a command line interpreter aimed to be very similar to MATLAB, but there is also a C++ API available for use. Since you refer to signal processing in contrast to image processing, I assume you mean audio processing, so you might need to look into the "signal" and "audio" packages in octave (which provide e.g. wavread). Keep in mind however, that the API is primarily designed for internal use in the octave CLI, and that it is considered not stable (i.e. the bindings might change). Also, since the library is not primarily meant to be linked with, the documentation is poor to say the least.

Alternatively, there is also GNU Scientific Library, which is not primarily aimed at signal processing, but provides complex numbers, matrices, fft and so on. On Windows, you would have to compile that yourself using Cygwin.


What exactly do we do in non-rigid registration?

I have found this survey paper very useful:
http://ralph.cs.cf.ac.uk/papers/Geometry/Registration.pdf

It gives a categorized overview of the field. I think section 3 and 4 are the most relevant ones for you.


finding frequency range and dominant frequency of a stored signal

To detect the fundamental frequency of your signals, you can use self-correlation. The maximum of the self-correlation will correspond to your fundamental frequency (the second maximum, as self-correlation is always maximum at 0). You can then interpret the lag value of this secondary maximum to find the periodicity of your signal, and the amplitude of the maximum to evaluate the strength of periodicity.
note : The value of the self-correlation at lag=0 is the energy of your signal. This can help as a scale for the periodicity strength.

Otherwise, to visually analyse the spectral content of your signal, you can just apply a Fourier Transform to it and plot the magnitude of the Fourier Transform.


Is it possible to combine two Sobel kernels into one?

A Sobel edge detector compares an approximation of the image gradient $|\nabla f_c(x,y)|$ to a threshold to decide if a pixel is an edge or not. A proper threshold must be determined and computed so that the comparison produces useful results. The approximation is given by the following:
$$|\nabla f_c(x,y)| \longrightarrow \sqrt{ (f \star h_x)^2 + (f \star h_y)^2 }$$
where $f = f[n_1,n_2]$ is the discrete-time image sequence, and $h_x = h_x[n_1,n_2]$ and $h_y = h_y[n_1,n_2]$ are the two directional Sobel kernels respectively.
As seen, the Sobel operator is nonlinear in the filters $h_x$ and $h_y$, and a single kernel, which would replace those two directional kernels, cannot be obtained in this case.


operations valid for predicting new codewords for a Linear block code having linear encoder mapping


MFCC, uncertain of my results and algorithm

For the computation of MFCC's you usually use a DCT-II. Just implement it yourself if you can't find a proper implementation anywhere. You only have relatively few log spectral values, and probably even fewer DCT coefficients that you will use, so I hope that computational complexity will be no big issue. This document is a nice overview of all the steps necessary for obtaining MFCC's.


Number of discrete samples required for the longest wavelength

If it is a real signal, as you are doing, then the negative spectrum will be the complex conjugate of the positive spectrum. Since there are N samples in frequency given N samples in time, and the frequency axis extends from 0 to N-1 where 0 represents DC and N represents the sampling rate (with the upper half of this from N/2 to N-1 equally representing the negative frequency axis) then yes in this case you need 2N samples to represent the spectrum without aliasing; and yes this is explained by the Nyquist-Shannon sampling theorem: https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem

However if you used an equivalent complex signal, then N samples would be sufficient.

$$x[n]=e^{j(30 \cdot 2\pi n/N)} + e^{j(60 \cdot 2\pi n/N)}  + e^{j(120 \cdot 2\pi n/N)} \\  + e^{j(260 \cdot 2\pi n/N)} + e^{j(330 \cdot 2\pi n/N)}$$


How is it possible plotting frequency vs. time in audio?

For traditional signal processing, time and frequency are dual variables, linked through Fourier transformations. This strong linkage yields a balance, sometimes called Heisenberg–Pauli–Weyl uncertainty inequalities. In everyday words, one could not be arbitrarily precise in time location and frequency determination. One cannot determine the frequency of a single sample.

However, in audio, meaningful information is usually carried by chunks of time, or collections of samples. So very often, people perform Fourier-like transformations over those consecutive or overlapping time chunks. They obtain information on how frequency contents evolve over time. A lot of specific methods, termed time-frequency analysis, have thus been developed, to provide  us with useful algorithmic tools to perform such analyses. They are legion: oversampled filter bank, spectrogram, Wigner-Ville, Choi-Williams, reassigned versions, S-transform, wavelets, etc.

In the tutorial: Time-Frequency Toolbox For Use with MATLAB (150 pages, no less), there are many comparisons. On a piecewise regular signal (sine then tone and linear chirp), they provide the following comparative plots that you can reproduce for your purpose:

One uses them  in audio for several reasons: to  analyze stationary parts, detect the onset of certains frequencies or harmonics, detect hearable information. Such techniques are at the very heart of mp3 audio compression for instance.


Meaning of lines on x- and y-axis of a centered FFT-spectrum of an image

For an image containing an oriented 1D sine pattern looking like a corrugated sheet like the ones displayed, the 2D Fourier transform magnitude spectrum generally contains information along a line (aligned dots). One indeed may expect (at least) three dots: one in the center for the average or DC component (which could be zero if the mean is subtracted from the image), and two symmetric ones whose distance is in inverse proportionality with the frequency of the sine.
Since the image is limited in size (hence windowed) and the spectrum is computed using an FFT (which hypothetizes periodicity of the input), one may observe "ripples" (corresponding with a cardinal sine kernel in the frequency domain, dual to the 2D rectangular window) and border effects (with a centered cross).
You can find additional descriptions in Fourier Transform or
slides 8 to 14 of 2D Fourier transforms.


A matlab implementation is in this answer: https://dsp.stackexchange.com/a/14201/5737

1) The wikipedia formula is a little bit too general.

2) If you know the basics of the with Fourier transform then you empirically know that: The image is viewed as being formed by superimposing a series of sinusoidal waves of various frequencies oriented in all kind of directions. Each "pixel" in the transform tells us the "intensity" of such a wave. The position of the "pixel" tells us the frequency and orientation of the wave. In practice, one wants to select only certain waves, having a specific frequency and a specific orientation.

So there you have it: The Gabor transform is one of many so called band pass filters that allows you to "cut" the Fourier transform and isolate only specific information. Another important information is that each Fourier "pixel" is a complex value (real and imaginary part)

3) Parameters: Two parameters are already shown:

3.a) The tuning frequency $f_0$, or tuning period $P_0$, or $\lambda$ establish to what kind of sinus wave the filter will respond best. ( $f = 1/P_0 = 1/\lambda$ or $f=\pi/\lambda$ depending on the specific implementation ) Basically, a smaller $P_0$ means a denser sinus wave. A larger $P_0$ means larger waves. $P_0$  is in pixels (3, 5, 30, etc. pixels) Don't go under 3 pixels or beyond $W/2$ or you get nasty effects. $W$ is the width of the image, if the image is square. You specify this using P0 parameter in the matlab code.

3.b) Central angle. These waves can have any directions. You want to select only waves at specific angle. So, the second parameter is the tuning angle, $\theta_0$ or $\theta$ in your formula. Usually, in radians. Orient in matlab code.

One cannot isolate only a certain frequency or a certain orientation. (Search for uncertainty principle in the textbooks. Yes, is similar to the one in physics) But one can tune how much of the nearby frequencies will leak. The next two parameters specify that:

3.c) $\Delta F$, frequency bandwidth expressed in octaves. Usefull values, 1.5, 2, 3. Larger values means capturing a broader range of frequencies. There is a price for a tighter bandpass, a poorer spatial localization. Why, again, textbook. FBW in matlab code.

3.d) $\Delta \theta$ the angle bandwidth. Expressed in radians. $\pi/3$ or $\pi/2$ works just fine. ABW in matlab code.

The relation between $\Delta f$, $\Delta \theta$ from matlab code and $\sigma, \gamma$ from Wikipedia have a formula but is not essential for understanding Gabor. $\psi$ is again not important for basic understanding.

So you have it: $X$ from Wikipedia is a 2D matrix of numbers containing a convolution mask. You take the original image, filter it with the convolution mask and get another image. This new image is the "Gabor response" for the original image.

The matlab code does these two steps together. Constructs the Gabor filter with the specified parameters and performs convolution. The results  ReConv,ImConv are the responses. Each  response "pixel" have a real and an imaginary part. If you want to use this code you usually must compute the energy for the response: $E = \sqrt{a^2+b^2}$ where $a$ is the real part of the response (ReConv) and $b$ is the imaginary part (ImConv), for each pixel.

There you have it:
1) Build a Gabor filter specifying $P_0, \theta_0, \Delta f, \Delta \theta$
2) Convolute your image with the filter. You will get two values per each pixel.
3) Compute the energy $E$ and get the intensity of the response for each pixel in the original image

Another intuition: Suppose you want to select edges stretching on an orientation perpendicular to $\pi/6$ and a certain width of 20 pixels. You can build a Gabor filter with $P_0 \approx 20, \theta_0 = \pi/6, \Delta f = 2, \Delta \theta = pi/2$. TAKE CARE that there is no equal relation between the dimension of your edge and the $P_0$ parameter. You should try various values and see what works best. The other parameters ($\Delta f, \Delta \theta$) touch them when you have some experience tuning the first two.

Hope it helps!

Cristi

Update

Here is a site which allows you to play a bit with Gabor parameters and note the results: http://www.cogsci.nl/pages/gabor-generator.php

After a quick look, both freqnency and angle bandwidths are tied together as "Standard Deviation in pixels ... to the Gauss envelope". The rest of the parameters are easily identified. However, note that there are numerical differences! Eg. the phase can be expressed in 0-1 interval or 0-$\pi$ interval.


relationship between block size and detected frequencies in goertzel algorithm

Since Görtzel computes a part of the DFT, the relation between detectable frequencies is the same. So with 8000Hz sampling rate and 600 samples you cover the frequency range of -4000Hz to 4000Hz with 600 frequency intervals, and both your frequencies are points of the subdivision, thus both should be perfectly detectable.

Can you show the code for the Görtzel iterations?

Or even before that, did you build the test signals using $\sin(\pi\cdot f\cdot t)$ or $\sin(2\pi\cdot f\cdot t)$? Only the second gives the correct frequency, the first could explain your observations, since the 600 samples would then contain $12.5$ periods of the oscillation for $f=333.33Hz$, but full $5$ for $133.33Hz$.


Is possible to perform algebric operation with bitwise operations in PDM domain?

Yes, this is sort of possible.   While PDM signals are full of high-frequency noise, at lower frequencies the "density" domain is limited to signals representing values from -1.0 thru +1.0.   This makes it difficult to implement true addition/subtraction since 1+1 exceeds the range of possible outputs, however if the output signal is scaled back into the allowable range, so averages like 1/2(A+B) and 1/2(A-B) can be calculated by multiplexing the two bitstreams into a single output stream (at twice the bitrate).

The best explanation of this I've seen was in a series of five articles titled "Signal Processing in the Density Domain" published in 2011 by Dave Van Ess in Electronic Design magazine.   The web copies of these articles exist but the equations and figures from the original versions seem corrupted by the web's evolution over time.   I'll include this link to the finale, Part V which wraps up how addition/subtraction and multiplication work.

While not easy, the missing figures can be retrieved from the internet wayback machine.


Why the following discrete type system,variable accumulator, is time invariant?

You are overloading the symbol $n$ which doesn't mean the same everywhere in your calculations.  Instead, write down what the sequence of input values is and the corresponding sequence of output values when the input is $x$.  DO NOT use the symbol $n$ anywhere ($n_0$ is OK to use).  So, your answer should look like

\begin{align}
y[-2] &= \cdots\\
y[-1] &= \cdots\\
y[0] &= x[-n_0] + x[-n_0+1] + \cdots + x[0] + x[1] + \cdots x[n_0]\\
y[1] &= x[-n_0+1] + \cdots + x[0] + x[1] + \cdots x[n_0+1]
\end{align}
etc.

Now, shift $x$ by $k$ places (say $k=1$ for starters) and call the result $\hat x$. Repeat the above computation for the input $\hat x$ and call the result $\hat y$.  Is $\hat y$ the same as $y$ except shifted by $k$ places? Repeat for $k+1$ etc until you have managed to convince yourself that regardless of what $k$ you choose, shifting $x$ by $k$ places results in $y$ also shifting by $k$ places.


When $\log_2(M)$ is even, the relationship betweeen BER and QAM if an AWGN channel is given as below (Assuming Gray coding so that a symbol error in most cases is one bit error):

\begin{align}
k&=\log_2(M)\\
y &= 10^{\frac{SNR}{10}}\\
P_e&=\frac{4}{k}\frac{\sqrt{M}-1}{2\sqrt{M}}\textrm{erfc}\left(\frac{\sqrt{3k\frac{y}{M-1}}}{\sqrt{2}}\right)
\end{align}

Using the above formula to solve for $y$ given your SNR and constellation $M$ should achieve the result you are looking for.

This is an excellent explanation for the derivation of the formula, specific to QAM.

And here is a nice summary by the same author for BER vs SNR of many different modulation schemes


How to convert the mean and variance of a processed received signal into a SNR or BER?

I also feel that the question needs a little bit more clarity but with what is stated so far, I think I can have a go at a response and edit it later if more information is added to the question.

  What I'd like to do is...

I do not intend to 'put you on the spot' with my question here: Have you checked how other people have tackled this? It seems that what you are trying to do is get an estimate of received signal quality (?). One of the simplest ways to achieve this is Received Signal Strength Indication (RSSI).

The estimation of both Signal To Noise Ratio (SNR) and Bit Error Ratio (BER) imply the use of a 'reference signal'. We will come back to this.

  What's going on is when we get our signal, we multiply it with the previous bit to get a correlation, and from that we sum it with the correlation from the previous bit to get our detection signal. This detection signal is what I'm taking the mean and variance of. If the mid-bit point of this detection signal is above a threshold (0 for now) then we say the bit is a 1, if it's below the threshold, we say it's a 0. I'm taking the mean and variance of all the 1's and the mean and variance of all the 0's, and am trying to come up with an average SNR type quantity.

This sounds like a differential codec (e.g. this one). Two questions emerge naturally here: The first one is, how do you know that the bit your are comparing against was decoded correctly, in order to re-use it to guide the decoding of your next bit (if I am getting the scheme correctly here) and the second question takes us back to RSSI and is "If you are trying to estimate statistics of the received signal timeseries with the objective of estimating the SNR, then why don't you use the RSSI?".

OK, so, let's see how are we going to do this. What is Signal to Noise Ratio?. It is the ratio between the strength of a "desired" signal versus the strength of an "undesired" signal (obviously, in the same context!). What is our "desired" signal here? It is the output of the modulator which is coupled to our medium somehow (antenna, speaker, light, whatever). What is the "undesired" signal here? It is the level of background "noise", undesired perturbations depending on the medium, which distort reception.

So, a way of guessing how bad our reception is going to be would be to "listen" to the background noise when there is NO reception (for some time) and create the N and then compare this with the strength of the signal S when we know that the transmitter is transmitting. Here, we form a ratio between a reference (our background noise) and a signal (received signal strength). The closer these two numbers are, the harder it is to accurately demodulate and eventually decode the signal.

If you absolutely must derive an SNR and a BER estimate, then you can examine the use of pilot signals or sequence. Communication systems that use a carrier signal, usually have a Phase Locked Loop (PLL) to lock on the carrier and track it. To allow for the PLL to achieve "lock", prior to transmitting the actual "information" signal, you can transmit a "dummy" sequence of some length, for the PLL to lock (e.g. listen to this clip). After the PLL "locks", the local oscillator (of the demodulator) is synchronised with the remote oscillator (of the modulator, at the transmitter side) and we can now start demodulating symbols. BUT we now do not know, in this bitstream, where does the actual information starts and where does it end. The pilot sequence may have been 100 symbols long. Maybe the PLL locked after 4 cycles, in which case, we have to throw away 96 "dummy" symbols, but maybe it locked after 50 symbols, in which case we have to throw away 50 "dummy" symbols. The point is that we don't know this 4 or 50. All we know is a line coming out of the PLL that is raised to tell us "I have now achieved lock".

So what we do, after the pilot tone, is to insert a start sequence. For example 01101001100101100. Once the system detects that, it knows that it can start decoding the actual information.

So, what you can do, is insert a long pilot sequence the structure of which you know (reference signal). For example 10101010101010101010101010101010101010101010101010101[BLOCK_START_SEQUENCE]. Once your PLL achieves lock, you don't throw away the symbols but you store them until you find the start sequence at which point you stop. Now, the pilot sequence could be 100 symbols long, the PLL achieved lock after 20 which means that we have decoded 80 symbols until the block start sequence. You know that these 80 symbols should be 10101010101010101010 (reference signal) and you count how many of them have been decoded in the right way. For example, 80 symbols decoded, 40 of them decoded right, that is a bit error rate of 0.5. For more information, you might want to see this paper.

Presumably, you can then use the BER in some adaptive scheme to select a strong or weak code. In addition, you might also want to look at equalisation and/or adaptive coding

Hope this helps.


Fragmented recovered image using a Phase Retrieval Algorithm named Error-Reduction

I think the effect you're seeing can be explained in terms of the shift theorem of Fourier transforms and the periodicity of the discrete Fourier transform.

It looks like the recovered object has been shifted upwards, and because the DFT is periodic, the top part of the object that's left the field of view has wrapped around to re-appear at the bottom of the field.
The reason the object is shifted is hard to see without knowing more about how your implementation of phase retrieval.  I can say that, due to the shift theorem, we know a shift in the object domain causes a linear phase in the Fourier domain, but does not affect the Fourier modulus.  So, the Fourier modulus alone can't help you find the correct shift.
I suspect you're leaving out some key information when you say "The only information given is the Fourier Modulus", because if that were strictly true you wouldn't be able to do phase retrieval at all.  As discussed in the paper, you need to have some other piece of information, typically the modulus of the object or the object support (the object support is a closed boundary that you know contains the object).  If you know the modulus of the object or at least its support, then reinforcing that constraint on every iteration of the error-reduction algorithm would prevent the object estimate from shifting.

As to your third question, I would point out that error-reduction is only the simplest algorithm discussed in the paper.  Gradient searches or the input-output methods described in the later sections are usually faster and/or more reliable.  In more recent years, successful methods that leverage nonlinear optimization algorithms such as BFGS have also become common. If you're interested, you might check out Prof. Fienup's list of publications.


Amplitude and Phase spectrum of a signal

When the spectrum is complex conjugate symmetric around $f_s/2$ (and multiples etc.) the corresponding signal is real valued and vice versa.

This holds for your example, meaning your signal should be real-valued. Your phase spectrum is interesting, since zero phase shift in all frequencies means your spectrum is real valued as well.


Understanding the transfer function of an FIR filter

Looks like a notation or transcription error. The "z" that's used as the transform variable to calculate the Z-Transform can not show up as an independent variable in the time domain as well. One of them needs a different name.


According to https://en.m.wikipedia.org/wiki/Decibel_watt

dBW is relative to 1 W and dBm is relative to 1 mW.  I’d just make sure to provide the correct units, but maybe you could ask your teacher if it matters.


Is there a method to calculate the maximum gain of a biquad?

in this previous answer, the squared-magnitude of a biquad frequency response is shown to be:

$$\begin{align} f(\phi) \triangleq |H(e^{j\omega})|^2
& = \left| \frac{b_0 + b_1 e^{-j\omega} + b_2 e^{-j2\omega}}{a_0 + a_1 e^{-j\omega} + a_2 e^{-j2\omega}} \right|^2 \\
& = \frac{B_0 + B_1 \phi + B_2 \phi^2}{A_0 + A_1 \phi + A_2 \phi^2} \\ \end{align}$$

where $ \phi \triangleq \sin^2\left(\frac{\omega}{2} \right) $

and

$\begin{align}
B_0 & = \left(\frac{b_0+b_1+b_2}{2}\right)^2 \\
B_1 & = - (4b_0b_2 + b_1(b_0+b_2)) \\
B_2 & = 4b_0b_2 \\
A_0 & = \left(\frac{a_0+a_1+a_2}{2}\right)^2 \\
A_1 & = - (4a_0a_2 + a_1(a_0+a_2)) \\
A_2 & = 4a_0a_2 \\
\end{align}$

maxima and minima occur when the derivative is zero.

the derivative of the above is zero when

$$ \phi = \frac{A_2B_0 - A_0B_2 \pm \sqrt{A_0^2B_2^2 + A_2^2B_0^2 - 2A_0A_2B_0B_2 - A_0A_1B_1B_2 \\ + A_0A_2B_1^2 + A_1^2B_0B_2 - A_1A_2B_0B_1}}{A_1B_2 - A_2B_1} $$

then evaluating $f(\phi)$ at those two values of $\phi$ is another laborious pain in ass.  and it won't be my ass.


Why only take 205 sample when using a sample rate of 10khz or 8khz?

The frequency of DFT bin centers is i*Fs/N, for integer i.  Thus different N will create different arrays of frequencies.  For some N, this will produce an array that has the lowest nearest bin center distance from the frequency set that you wish to measure (in thus case, a vector of DFT frequencies).  According to the referenced document, for some standard sample rate, N=205 will produce that vector with the closest DFT bin centers to DFT frequencies.

Since each DFT bin has a Sinc shaped roll-off, the closer a bin frequency is to a desired measurement frequency, the lower the magnitude measurement error will be without needing further interpolation computation.

For N greater than 320, the DFT will have greater frequency resolution, but the time resolution gets worse and will overlap multiple tones, which creates greater ambiguity in which tone is being detected.


Is there an a method to fit a wave created from two wave?

As a precautionary statement, the FFT will give you coefficients of the frequency components closest to the integer multiples of $1/T$ in Hz where $T$ is the length of the sequence in seconds being used for the FFT. If the components are not on these integer boundaries, then you will get several results (called spectral leakage) where the actual components (if information is known such as in the OP’s case of having just two tones) can be estimated as being in between two large adjacent results.
Further when combining the waveform you must pay attention to the phase as well as amplitude of the result in order for it to match the time domain waveform shown. The FFT will provide both a scaled magnitude and phase. The FFT results are the coefficients for complex frequencies (of the form $e^{j\omega t}$ but if it is known that the waveform is real then you can use the first half of the FFT as representing the coefficients of real frequencies in the form of $\cos(\omega t)$, scaled by $N/2$ where $N$ is the number of samples in the sequence used for the FFT.


Relation between Wireless Power Transfer, Wireless Information Transfer and the definition of sum rate

"sum-rate" is just the sum of all rates of communications happening in a network.

The statement is pretty self-evident: The more time your nodes can communicate, the more data they can transfer.


If with "Comb Function" you are referring to the periodicity of the spectrum of a sampled signal, first see Intuition for sidelobes in FFT for further details on the unfolded view of a digital spectrum (specifically how it repeats over the range of 0 to $F_s$ where $F_s$ is the sampling rate.)

(If you are referring to the use of "Comb Filters" then see this explanation of the CIC Decimator: Fast Integer 8 Hz 2nd Order LP for Microcontroller, and I can elaborate upon request on the CIC Interpolator which has a similar structure).

Upsampling and interpolation in general can be viewed as diagrammed in the figure below (once understanding the unfolded view of the digital spectrum from the top link above).

The frequency is upsampled by a factor of M by inserting M-1 zeros. In the case of the example in the figure, 3 zeros are inserted between every digital sample to increase the sampling rate by a factor of 4.

What happens when we do the zero insert is we get our original spectrum intact, along with new frequency components (replicas of our original spectrum) that now exist in our primary digital frequency span of $0$ to $F_s$. These are the exact same components that existed adjacent to all multiples of the sampling rate as described above in the unfolded view of the digital spectrum. If we can create a digital filter (Interpolation FIR or IFIR) that will pass the signal of interest with no distortion, and completely eliminate the unwanted components that would otherwise be distortion, as shown in the figure below, we will have a perfect interpolation between the original samples (the filter will grow the zero values to the perfect interpolation value needed, by nature of having all distortion frequencies eliminated while not distorting our spectrum of interest). Of course such an ideal filter is not realizable, but we can control and decide on how good of an interpolation we want based on the complexity (number of taps) in our filter design.

Also one may ask why the zero insert, why not do a first order hold (meaning use the current value as the value to insert instead of inserting zeros)?

To demonstrate the considerations with this, consider that a zero order hold is done mathematically in the time domain by convolving our waveform with a digital "pulse", while a zero insert is done by convolving our waveform with a digital "impulse".

Convolution in the time domain is multiplication in the frequency domain, so to explain using analog terms which the digital system will approximate- for the case of a first order hold, we are multiplying our desired frequency response by a sinc function (FT of a pulse) leading to passband droop, but for the zero insert, we are multiplying our desired frequency response by a constant (FT of an impulse), leaving our original passband undistorted.   The plots below demonstrate this for our intepolate by 4 example with the exact digital solution of what would occur.

Convolving with the digital equivalent of a sinc function is not a terrible way to go, and in fact its nulls will be right in the center of the frequency bands we want to suppress (very convenient). It is just important to understand the implications of passband droop and how to compensate for that if you take that approach (see how to make CIC compensation filter for a very simple 3 tap "inverse Sinc" compensation approach I have used). The CIC Interpolator is in fact such a "first-order hold" interpolation method that is very popular because of its simplicity.  Otherwise is absolute minimum distortion of your original waveform over its full usable bandwidth is important then zero-insert with a good filter design based on defined passband and stop band requirements is the way to go. Filter design is optimized by only rejecting the bands where the spectrum is replicated. Matlab's firls and firpm filter commands (Least Squared and Parks-McCLellan) both support creating multi-band filters which are perfect for this purpose. Further, if you go down the path of zero insert and optimum filter design, the end result can be implemented efficiently as a polyphase interpolator as shown in the figure below.  The process is to design the filter as described above, and then map the same coefficients into the polyphase structure as shown in the figure. The figure is descriptive of the mapping process in showing only 8 taps, but of course with only 8 taps in the original filter you cannot expect very good performance. An actual interpolate by 4 filter may have 64 taps for example, leading to 4 polyphase filters each with 16 taps, mapped the same way from the original base filter to the polyphase structure. The more taps the better the base filter can minimize pass band ripple and maximize the stop band rejections.

The reason for this mapping is very clear by following what occurs with the zeros that were inserted in the original filter- they mask several coefficients that are not needed on any particular given cycle (see figure below):


It's instructive to see what an ideal filter adding a constant phase shift would look like. If $\theta$ is the desired phase shift, the corresponding ideal frequency response is

$$H(e^{j\omega})=\begin{cases}e^{-j\theta}&,\quad 0<\omega<\pi\\
e^{j\theta}&,\quad-\pi<\omega<0\end{cases}\tag{1}$$

Using the sign function $\text{sign}(\omega)$, the frequency response $(1)$ can be rewritten as

$$H(e^{j\omega})=\cos\theta-j\,\text{sign}(\omega)\sin\theta,\quad -\pi<\omega <\pi\tag{2}$$

With the DTFT correspondence

$$-j\,\text{sign}(\omega)\Longleftrightarrow g[n]=\begin{cases}\frac{2}{\pi n},&\quad n\text{ odd}\\
0,&\quad n\text{ even}\end{cases}\tag{3}$$

the impulse response corresponding to $H(e^{j\omega})$ is given by

$$h[n]=\cos\theta\cdot\delta[n]+\sin\theta\cdot g[n]\tag{4}$$

where $g[n]$ is the sequence on the right-hand side of Eq. $(3)$, which is the impulse response of an ideal discrete-time Hilbert transformer. Eq. $(4)$ shows that an ideal phase shifter can be implemented as a weighted parallel connection of a wire ($\cos\theta\cdot\delta[n]$) and a Hilbert transformer.

So your problem can be solved by using any of the many available designs of discrete-time Hilbert transformers. Note that you can get much better performance for a given filter order by taking into account that the approximation needs only be accurate in the given frequency band. For a frequency-domain design method this just means that in the formulation of the desired response given in $(1)$ you replace the positive frequencies by the frequency interval of interest and leave the rest as a "don't care" region. The same is done for the negative frequencies.


How to count the number of threads in a bolt or nut using image processing?

The most important thing in this problem is to be able to control the scene and more importantly:

The lighting
The camera position & angle

The Lighting:
Simply put, the lighting is what creates the essential contrast to delineate the thread. If you light the thread "head on" (light from the viewpoint of the camera), there will be no contrast or "ups and downs" to count. The better position for the light is at the top and slightly forward of the face of the thread so that a) All thread peaks are lit (both those that are on the top, but also those that are below them) and b) The top thread casts a bit of shadow to the one below it. Try to think of the thread as a difraction grating, only here, you need the shadows.

The Camera: The camera distance from the thread and angle to the thread are very important, especially if you wish at some point to derive the pitch of the thread. To derive the pitch of the thread, you would have to count lines per unit of distance.

The Camera Distance: The only way to relate the pixel size to a physical unit of length in a monocular image, is to fix the distance between the subject and the lens, similarly to how they do it in archaeology photography. There, we assume that the ruler's plain coincides with the focus plane of the lens, or is too near to be creating a negligible error.
The Camera Angle: The camera angle is incredibly important in this problem because it is responsible for an optical distortion that would affect your measurement. In the nut example, it is impossible to image the thread vertically ("head on"), because you are trying to image the inside of a very small and thick cylinder. Therefore, you have to look at it from an angle. In the bolt example, it is impossible to image the thread vertically ("head on"), because there will be no shadows to create contrast. A pen, imaged at an angle will appear shorter than it really is. Similarly, the thread imaged at an angle, will not give you equidistant peaks. Specifically, as the thread "screws away" from the camera, its pitch will be increasing because the lines that denote the threads would appear to get closer and closer and closer. Therefore, the angle that the camera is set against the screw determines a) The error of measurement and b) the limits of measurement (how fine / coarse a thread you can measure).

AFTER determining / fixing these parameters which are very important to the problem you are dealing with, you can start thinking about extracting the thread. If you get the imaging conditions right, the thread extraction via image processing will be much more easier and straightforward than trying to do it over any random image.

Trivial way of extracting the thread:

Assuming that the thread is aligned and centred within the camera's Field Of View (FOV), you are using something like Octave and the image has been loaded in I with something like I = imread('someNutImage.jpg');:

Threshold the image

Because the thread is (supposedly) aligned with the FOV, take the middle column of values. That would look something like:

Notice here how the pitch changes visibly because of the camera angle.

Now, "count" the gaps. This is a trivial "pulse counting" excercise which I have outlined previously. Prior to counting the pulse widths, you need an absolutely saturated signal. Therefore Q = double(I(:,320)>128).*255; (Where 320 is some column image, I am taking the middle here as explained above). Now, differentiate the signal to get the "pulse fronts" and then differentiate that again, to get the "inter-pulse-front" distances. That is:

pulse_fronts = find(abs(diff(Q)));
pulse_lengths = diff(pulse_fronts);

For the signal I am using here, I get:

1, 24, 64, 11, 8, 13, 6, 14, 6, 13, 6, 13, 6, 13, 5, 12, 6, 11, 5, 12, 5, 10, 6

That is, the signal starts with a spike pulse of length 1, followed by 24 zeroes followed by 11 ones, followed by 8 zeros, followed by 13 ones and so on.

Of course, you are not going to "count the gaps" manually. At this point, you can do a histogram and obtain the two most frequent values. For example, with a hist(pulseLength,max(pulseLength)), you get something like:

The "low" most frequent value (6) are the zeroes and the "high" most frequent value (13) are the ones. Once you have these two values (from the histogram), you know the full length of your thread in pixels. To convert it to mm, you need to determine the pixel size.

Less trivial way of extracting the thread:

The above method allows you to filter (if required) the pulse lengths to get rid of various artifacts. Another way to extract the thread would be to use the Discrete Fourier Transform of Q to assess the frequency of the pulse train.

For example: semilogy(abs(fft(Q))), gives you something like:

This image is 480 pixels tall. This means that the highest frequency (in terms of alternating black/white lines per unit of length) is 240. The first peak of the DFT occurs at 26 lines. From the method above we had ~13+~6=~19. The DFT will give you an average pitch over the whole acquired signal. Again, to give these numbers a physical meaning, you would need to know the pixel size which is only known by fixing the distance from the camera.

Even less trivial way of extracting the thread:

Once you threshold the image, either apply thining or a simple high-pass convolution mask to preserve the edges and then a Hough Transform. This will generate one "dot" for every slanted "line" that is created by the thread. What you would be looking for then would be a series of dots, one next to the other signifying the thread. By measuring the distance between these dots, you can extract your thread information. This is probably the most accurate method because the Hough Transform is very sensitive and accurate at recovering lines from very noisy images like the thresholded one we are talking about here. It is also the most computationally expensive option.

Each one of these methods has pros and cons for estimating the thread of a nut or bolt.

Finally, the most important part in such applications is to fix the imaging geometry so that you know what you measure and what sort of limits this measurement has.

Hope this helps.


Depends on where you're trying to implement/the resources available to you.

I've used Support Vector Machine learning for that task; you have to feed in individual characters to train your algo.

take a look at python's scikit-learn and if you ned to implement on FPGA check out: MyHDL


Graph Fourier transform: the adjoint notation for the eigenbasis matrix

You terminology is something I am not used to.

$$ \hat{x} = \mathcal{F}_{G} x $$

If by $\hat{x}$ you mean the DFT (complex except very special cases), $x$ is your input (real or complex), and $\mathcal{F}_{G}$ is the matrix composed of the sinusoidal basis vectors, then your assertion that $U$ is real is dubious.  $\mathcal{F}_{G}$ is most certainly not.  With the proper normalizattion:

$$ \mathcal{F}_{G}^{-1} = \mathcal{F}_{G}^{*} $$

Now, for a known $x$ where $\hat{x}$ is known to be complex, how can $U$ be real?

Since they are using the complex conjugate, I would say that U is expected to also be capable of being complex.  I have not read the paper so I don't know where you picked up that assumption.


A)

The quantization of the photocurrent is an actual physical phenomenon.

I'll use $\Delta T$ in place of your $dT$. It is not a variable describing the actual physical process but rather the time step of your approximate model. Effectively you are filtering the continuous-time Poisson point process by a continuous-time "moving average" filter with a rectangular impulse response of length $\Delta T$ and sampling the filtered signal uniformly with a sampling period $\Delta T$. The filter acts similarly to a low-pass filter, but is not very effective at band-limiting the process. Therefore you have aliasing and it becomes difficult to interpret your results. If you don't know what aliasing means, have a look at the question: "explanation of Hybrid systems?" and its answers for an introduction on sampling of continuous-time signals.

Sampling the RC-filtered Poisson Process

Non-uniform sampling: The impulse response of a first-order low-pass filter is a right-sided exponential function. It would be easy to convolve the Poisson point process directly by it. Non-uniformly sampling the filtered signal at the event times would give an unambiguous representation of the signal that allows exact reconstruction of the signal from the samples. This would not be sampling of a band-limited signal but would take advantage of the special form of the signal. See math later below.

Uniform sampling: If the cut-off frequency of the first-order low-pass filter is low enough relative to the sampling frequency, then uniform sampling of the filtered signal does not cause much aliasing and you will have quite a clean discrete-time signal (time series) to work with. When time stepping, the signal decays exponentially from the previous sample to the current sample, unless there have been events within the last time step. In that case you'd calculate the contribution of the filtered events and add that into the current signal value. That's the beauty of the right-sided exponential impulse response: it allows to forget when the past events occurred, combining their contribution to a single value that keeps decaying exponentially. The drawback of uniform sampling is that if there are two or more events between successive samples, then their exact times cannot in general be reconstructed from the samples of the filtered signal.

There are two easy ways to generate the event times of a Poisson process: 1) The time between successive events follows an exponential distribution, so you just need a pseudo-random number generator with an exponential distribution, and generate the times of the events in order until you are at the simulation end time. The time of an event is obtained adding the exponentially distributed random number to the time of the previous event. Note that it is only a coincidence that also the low-pass filter impulse response is a one-sided exponential function. 2) Generate a Poisson distributed random number to get the number of events over the simulation length, pick that many event times from an uniform distribution that is as long as the simulation and sort the event times. Or you can generate event times for each sampling interval individually.

In this presentation I will not use units, for my own convenience. The impulse response $h$ of a resistor-capacitor (RC) filter, as function of continuous time $t$, is of form:

$$h(t) = \begin{cases}0&\text{if }t<0,\\\displaystyle a e^{-at}&\text{if }t\ge0,\end{cases}\tag{1}$$

Figure 1. RC-filter impulse response (Eq. 1), with decay parameter $a = 1$.

where $a$ is a positive constant that determines the decay rate. If this is convolved with a realization of a Poisson process with $K$ Dirac impulse events at times $p_k$, $0 < k \le K$, then we get the signal $y$ of which you can find some literature under the keyword filtered Poisson process:

$$y(t) = \sum_k^{K-1}h\left(t-p_k\right)\tag{2}$$

Figure 2. A realization of an RC-filtered Poisson process (Eq. 2), with $a = 1$.

We can set up a sample point at time $t_n$. If we know the signal value  $y(t_n)$ at that point, then we can forget the events at or before that time and still reconstruct $y(t)$ at later times $t \ge t_n$:

$$y(t) = y(t_n)\,e^{-a(t-t_n)} + \sum_{k=0}^{K-1}\begin{cases}0&\text{if }p_k \le t_n,\\h(t-p_k)&\text{if }p_k > t_n,\end{cases}\quad t \ge t_n,\tag{3}$$

or in a less formal short-hand notation:

$$y(t) = y(t_n)\,e^{-a(t-t_n)} + \sum_{k,\,t_n < p_k}h(t-p_k), \quad t \ge t_n.\tag{4}$$

This is a property unique to the one-sided exponential function $h(t)$. The effect of an event at time $p_k$ is already included in $y(p_k)$ with the way Eq. 1 is currently defined. We can make explicit also that because the impulse response is causal, the events happening after time $t$ do not affect $y(t)$:

$$y(t) = y(t_n)\,e^{-a(t-t_n)} + \sum_{k,\,t_n<p_k\le t}h(t-p_k), \quad t \ge t_n.\tag{5}$$

You can do non-uniform sampling at the event times by setting $t_k = p_k$, or uniform sampling with sampling period $\Delta t$ by setting $t_n = t_0 + n\,\Delta t$. In your data structure you'd just store $t_n$ (no need in case of uniform sampling) and $y(t_n)$. Python example:

import numpy as np
import bisect
import matplotlib.pyplot as plt

sim_length = 10  # Simulation length
rate = 100  # Event rate per time unit
a = 1  # Filter decay constant

K = np.random.poisson(rate*sim_length)  # Number of events
N = K + 1  # Number of samples. Always have one sample at simulation start
ty = np.empty(N, dtype = [('time', float), ('value', float)])  # Sample times t and signal values y
ty['time'][0] = 0
ty['value'][0] = 0  # Start signal at value 0
ty['time'][1:] = np.random.uniform(0, sim_length, K)  # Event times
ty[1:].sort(order='time')  # Sort events by time
for k in range(K):  # Sample signal at event times
    ty['value'][1 + k] = ty['value'][k]*np.exp(-a*(ty['time'][1 + k] - ty['time'][k])) + a

def sample_y(ty, t):  # Sample signal at arbitrary time t
    k = bisect.bisect(ty['time'], t) - 1
    return ty['value'][k]*np.exp(-a*(t-ty['time'][k]))

# Create a time series of the RC-filtered Poisson process
t = np.arange(0, sim_length, 1/10000)  # Delta t = 1/10000
y = np.zeros(t.size)
for k in range(t.size):
    y[k] = sample_y(ty, t[k])

# Plot
plt.plot(t, y)
plt.xlabel('t')
plt.ylabel('y(t)')
plt.show()

Figure 3. RC-filtered Poisson process plotted by the above Python script.

Figure 4. RC-filtered Poisson process reparameterized with rate = 1 and a = 10.


Programming a parametric equalizer - preventing clipping

There's no universally acclaimed correct design of equalizers. Just like there are multiple car designs per brand, there are also multiple techniques for achieving desired effects.

It's true that your channel filters will interact with each other and change the effective resulting gain at each frequency, but mostly pronounced at the adjacent channel peak though.

The amount of interaction depends on your selected DC gains (which are unity for your application) and channel bandwidths and bandwidth gains.

There are recommended practices to prevent distortion levels. Headrooms are included so that peak gains will not create clipping. Or you can allow some...?

Robert Bristow Johnson is a member here and he's one of the designers of those (cookbook) parametric equalizers that you want to use in your application.


Given the horizontal and vertical component of the distance vector from the source to the sensor array be X and Y, then the angle between the sensor array and the source is given by:
$$ \theta = \tan^{-1} (Y/X) $$

If the linear distance between each sensor is $d$, each sensor should experience a delay:
$$ \delta t = d \sin(\theta)/c $$

Now if you can measure the delays (high precision measurement will lead to higher accuracy), and have the knowledge of $d$, you can find the angle $\theta$.

Warning: Practical errors and considerations are not taken into account in this approach.


What you describe is indeed aliasing. If there are frequency components above half the sampling frequency then they will be folded back to the interval $[0,f_s/2]$. This folding back is called aliasing, no matter whether these aliased components overlap with other parts of the signal spectrum or not.
Note that you can also use aliasing to your advantage if you're careful that no overlap occurs. An example is bandpass sampling, where you deliberately choose a sampling frequency that is much smaller than twice the highest signal frequency.


Need help in understanding minimum variance estimator and CRLB concept using an example

You may want an estimator with minimun variance because each time you calculate it, it is likely to be closer to its expectation, and in the case that the estimator is unbiased, closer to the true parameter. An estimator is inconsisten if increasing the number of samples of the dataset to calculate it doesn't reduces its variance.

The sample mean estimator of X:
$$
\bar X = \frac{1}{N}\sum_{i=0}^{N-1}x_i
$$
is unbiased and its variance is:
$$
\mathbb var(X)=\frac{\sigma^2}{N}
$$
It is clear that as N increases, the variance gets smaller, so it is consistent.

Regarding question 1, NO. For instance, one of the two estimators can have a higher CRLB always, for both high or low SNR. The criteria for choosing an estimator depends on the problem, the prior knowledge about it, the computation constraints and so on. Normally it is desired a minumn variance and a non-biased estimator, i.e. an efficient estimator. But, if your decision is only based on the CRLB, then choose that one with the smaller one.

Take this 'silly' example: for a dataset you decide that the estimator for the mean is always (arbitrarily) 3.5 even if you know nothing about the data and you do no calculations... It is clear that this estimator has variance 0, because it will be always 3.5, but it is strongly biased, it is randomly biased indeed... is it a good estimator? No, although it has 0 variance. The key point here is that yo can't have an unbiased estimator with a lower variance than its CRLB (if the CRLB exists). If the variance of your estimator is lower than the CRLB, it is biased.

If you have some experience with computer simulation, Matlab or similars, it is a good exercise running Montecarlo simulations and plotting the variance and the CRLB of each realisation versus the SNR for a fixed number of samples, or/and for a fixed SNR plotting the variance and the CRLB of the estimator versus different dataset lengths.


Steady-State Output from Transfer Function

You may know that one important property of linear time-invariant (LTI) systems is that the complex exponential $e^{j\omega_0}$ is an eigenfunction, and the corresponding eigenvalue is given by the system's frequency response evaluated at $\omega_0$.

So the response to an input signal $$x[n]=e^{jn\omega_0}$$ is given by $$y[n]=H(e^{j\omega_0})e^{jn\omega_0}$$

Now try to show that from this it follows that the response to a sinusoidal input signal $$x[n]=\sin(\omega_0n)$$ is given by $$y[n]=\big|H(e^{j\omega_0})\big|\sin\left[\omega_0n+\angle H(e^{j\omega_0})\right]$$ (if the system is real-valued). Now should you have all the information to compute the reconstructed output signal.


What is the time-integration property in the Fourier series analysis?

Let's consider the complex form of the Fourier series of a $T$-periodic function $x(t)$:

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

with $\omega_0=2\pi/T$. The integration property says that the Fourier series of

$$y(t)=\int x(t)dt\tag{2}$$

is given by

$$y(t)=\sum_{k=-\infty}^{\infty}b_ke^{jk\omega_0t}\tag{3}$$

with

$$b_k=\frac{a_k}{jk\omega_0},\quad k\neq 0\tag{4}$$

where we require that $a_0=0$.

Note that in $(2)$ I used the indefinite integral, and I chose the integration constant to be zero, such that $b_0=0$. In some texts you may find the definite integral instead of the indefinite integral:

$$y(t)=\int_{-\infty}^{t} x(t)dt\tag{5}$$

The problem with $(5)$ is that for certain functions $x(t)$ the definite integral doesn't exist with $-\infty$ as the lower bound. Note that the lower bound in $(5)$ corresponds to the integration constant of the indefinite integral $(2)$. This becomes obvious if we consider the relation between the indefinite and the definite integral:

$$\int x(t)dt=X(t)+C\tag{6}$$

$$\int_{t_0}^tx(\tau)d\tau=X(t)-X(t_0)\tag{7}$$


Brief introduction to wireless communication and communication terminologies

Well, these are very fundamental questions, and I think they are written in any good book. But maybe it's hard to extract the exact information from them, which you ask for. So, let me briefly answer:

In wireless, the channel is usually modelled as a FIR filter with impulse response $h(\tau)$. But, this channel might change over time (due to e.g. mobility), to make it a time-variant impulse response $h(t,\tau)$. Moreover, the channel is often modeled as a tapped delay line, $h(t,\tau)=\sum_ia_i(t)\delta(\tau-\tau_i(t))$ where $a_i(t)$ are time-varying gains of the taps and $\tau_i(t)$ are the time-varying delays of the taps. The fact that $a_i(t)$ has complex values does not make it Rayleigh fading. Rayleigh fading describes the temporal distribution of $a_i(t)$. If over a long time range the distribution of each path gain is a Gaussian, it becomes Rayleigh fading. The doppler spread of the Rayleigh fading describes how quickly the value $a_i(t)$ changes over time. In the special case of block-fading $a_i(t)$ is a constant (the channel remains static) for the whole transmission. Then, when a new transmission starts, a new channel is generated. Then, the realizations $a_i$ follow a Gaussian, to make it block Rayleigh fading.

A multipath channel corresponds to the tapped-delay line model. Have a look at one of my articles about multipath propagation, maybe it helps. The number of inputs and outputs depends on the number of TX and RX antennas of your system (talking about MIMO systems here) or the number of users. For single-antenna systems there is just one input and one output.

The input to the physical channel is analog, if you model it as a tapped delay lines. However, you can assume transmission in passband, then all signals are real. Or, you can do an equivalent analysis in baseband, where all signals are complex but the math is easier to treat. For the description of transition from baseband to passband, you can look at another article of mine.
Though, there are other channel models such as BSC or BEC, though these are not used to model the transmission of an analog signal over the wireless physical channel, but abstract away from that.

Modulation is performed after encoding ;-). Modulation means to generate a signal which is adapted to the channel, a a function of the transmit bits. It can e.g. include mapping the bits to complex QAM constellation points.

Encoding takes your bits you want to transmit (the payload) and adds redundancy via a channel code. This way, if at the receiver side, one bit is lost, it can be recovered by appropriate decoding.


Why does the solution for the voltage of this RLC circuit result in a maximum capacitor value?

Why the fear of complex numbers? Just do a distinction of cases and generate solutions for the real and the complex case. $b-a$ in the denominator will always be real, as the imaginary parts will cancel out each other, so complex numbers will only be in the exponents, wich will lead to some nice sinusoidal functions, which is expected, as this is an oscillating system. The real case just represents the overdamped case with no oscillation. If you don't know what I'm talking about have a look at this wikipedia article:
https://en.wikipedia.org/wiki/Damping_ratio


What is the relation between the domain of Fractional Fourier transform and time and frequency?

It think it is more correct to say that the $u$ axis is a linear combination of the $x$ and $\xi$ axes:

$$(x, \xi) = (1, 0)x + (0, 1)\xi = \left((1, 0)\cos(\alpha) + (0, 1)\sin(\alpha)\right)u,$$

that is:

$$ x = u\cos(\alpha)\\
\xi = u\sin(\alpha).$$

The definition of the fractional Fourier transform from that Wikipedia page is:

  $$\mathcal{F}_\alpha[f](u) =  \sqrt{1-i\cot(\alpha)} e^{i \pi
\cot(\alpha) u^2}  \int_{-\infty}^\infty  e^{-i2\pi (\csc(\alpha) u x
- \frac{\cot(\alpha)}{2} x^2)} f(x)\, \mathrm{d}x.$$

If you use another variable, $y$, for the integral, and substitute $u = \sqrt{x^2 + \xi^2}$ and $\alpha = \text{atan}(\xi, x)$, then the fractional Fourier transform can be written as:

$$\mathcal{F}[f](x, \xi) = e^{i \pi \frac{\xi}{x} (x^2+\xi^2)} \sqrt{1-\frac{i \xi}{x}}\int^\infty_{-\infty}e^{i \pi \frac{y}{x} \left(\xi y-2 x^{\frac{1}{\xi}} \sqrt{\xi^2+x^2}\right)} f(y)\, \text{d}y.$$


I believe you are talking about soundwaves only, so I'll approach the answer from there.

The instant values a signal takes don't tell much about the sound it may generate. If a signal is zero at a certain instant $t_0$, it doesn't mean that at that moment there is silence. A single sample can't represent any sound because only oscillations can generate sound. I hope you see why a single sample can't possibly represent  an oscillation of any kind.

Take for example a perfect sinewave. The signal $x(t)=\sin(t)$ is $0$ for infinite values of $t$, but the sound you hear is nonetheless uniform and constant. Also, it takes both positive and negative values, but the note you hear is always the same.

If a signal is $0$ for a given period of time though, then for that interval there will be silence indeed. In fact, if a signal is constant for a given period of time (not necessarily $0$), you will hear no sound at all (but you would be wasting power nonsensically).

In a nutshell, the instant values of a signal represent no sound whatsoever. It's the change over time of these values that can represent and generate sounds.


Abstracting this a bit:

This doesn't look too bad for a signal!

So, first of all, you might be driving your sensor and digitizer too hard – artifacts due to nonlinearity would be at multiples of the fundamental tone! Idea is the following: If the transfer function of my system isn't $y = \alpha x$ but (also contains) a function $y=\alpha x^2$, then a sine in $x$ would lead to two new signals: one at 0 Hz, and one at twice the original sine's frequency. That's intermodulation as understandable by the classical trigonometric addition / multiplication theorems!

So, first thing: make sure your system is nice and linear over the whole range; for systems that take vibrations, exciting the thing with a large acceleration of known amplitude might be hard, so a first start would be plotting a histogram of the amplitude of your signal. Does it look nice and somewhat smooth with few values at the extreme values and most at low amplitudes?

If that's not the case, you might simply need to use a sensor / sensing system that can work with a larger range of values. That might just mean your sensor is working fine and linearly, but your ADC needs to get the values attenuated to map them to the whole range (you'd see that as sharp clipping if you plot your signal over time!); but based on your relatively benign "ghost signals", my guess is that your acceleration sensor simply isn't that linear. Maybe it's not meant for that range of acceleration, maybe you're not supplying it with a supply current or voltage high enough to represent these values? Maybe it's not even specified to do anything linear in frequency ranges above a certain cutoff frequency?

Another thing: Your plot goes up to what I presume is half the sampling rate, and I see no cutoff of any filter that was applied to the signal.

Are you properly anti-aliasing filtering your signal prior to sampling it? If your signal contains e.g. the oscillation of a MEMS fork at 5 kHz (wild guess!), but you're sampling at 4.8 kHz, then that oscillation would appear at 0.2 kHz, and a lot of multiples of 0.2 kHz. Do not mess with Nyquist unless you want aliases! Quite possible, the accelerometer actually generates a signal that contains the intermodulation of MEMS oscillation and signal to be observed:

$$\sin(2\pi f_{machine} t+\varphi_{machine}) \cdot \sin(2\pi f_{oscillator} t)$$

Which in turn would be

$$\frac12 \left(\cos(2\pi(f_{machine}-f_{oscillator})+\varphi_{machine})-\cos(\cos(2\pi(f_{machine}-f_{oscillator})\right)\text.$$

Staying with the number from my example, $f_{oscillator}=5\,\text{kHz}$ and your signal at 100 Hz, we'd thus should see components at 4.9 kHz and 5.1 kHz – except we don't, because we're not sampling fast enough, so at our 4.8 kHz sampling rate, these would end up at 100 and 300 Hz; and then, combine that with the non-linearity intermodulation, you get the funkiest frequency content.

So, first of all, if you know you're not expecting oscillations above 100 Hz, use an analog filter to get rid of anything above. If you still get signal components above that frequency, your filter isn't sufficient or your ADC is driven into nonlinearity.

On the other hand, without any education on the mechanical side of things: doesn't seem so unlikely to me that materials aren't perfectly rigid, and that something oscillating at 100 Hz might be "fluttering" at higher frequencies, especially if things like spring elements aren't the perfect linear springs I was told to assume in school physics.

In conclusion:

You're looking at something non-linear in nature. Whether it's your measurement system or the system you're observing is something we can't tell you, but systematic elimination of error sources will probably yield results pretty quickly.

Start with reducing sensitivity of your measurement system.

Also, make double sure you're properly anti-alias-filtering your signal prior to digitization.


Avoiding edge repeats in a FFT convolution

Zero-pad the edges of the image, then truncate the excess edge from the result


How to identify and differentiate frequency and time in EEG data using python?

It is hard to answer your question, since you do not seem to have experience with EEG data and/or general signal processing.

You say, your data is sampled at 200 Hz, which seems good to me for EEG data.

What does it mean by sampling in EEG data?
That means, that the (continous) EEG signal is choped up into discrete values. This is always necessary for further signal processing with computers.
I am getting the sample rate and the frequency when I plot the data as shown below. Shouldn't it be the time and frequency?

I don't know what you mean by I am getting the sample rate and the frequency when I plot the data. You certainly don't. The sampling rate/frequency is a fixed number given by your experimental hardware. You don't get it from your data - you have to know it from your experiment to get meaningful results from your data.

I also don't know what Shouldn't it be the time and frequency? means. Your EEG data is certainly a voltage signal over time, hence you have a signal in time domain.

Your sampling rate is given from your experimental settings and is 200 Hz as you say. Hence, each time point is separated by $$\Delta t = \frac{1}{f} = \frac{1}{200\,\textrm{Hz}} = 5\,\textrm{ms}$$

Your total signal that you have plotted is therefore $$4097 \times 5\times 10^{-3}\textrm{s} = 20.485\,\textrm{s}$$
long.

That is what you know from your data. Any further data extraction is up to you or your research question. If you'd like to get the frequencies that are in the signal, try to Fourier-transform it and extract the frequency region that you are interested in. As I see it, you have some apparent frequency components at approx. 13.4 Hz, 57.7 Hz, 81.0 Hz, and 88.0 Hz. Probably you are only interested in the lower spectrum around 13.4 Hz...


If you want to have the best SNR possible, you put a null on the signal that is leaking it's energy.  In order to do that, you either have to know the location of the interferer or use an adaptive filtering scheme.

If you want to stick with a fixed window, the Taylor window isn't a bad choice,  It's like a Chebyshev and rolls off further out.

A. Villeneuve, "Taylor patterns for discrete arrays," in IEEE Transactions on Antennas and Propagation, vol. 32, no. 10, pp. 1089-1093, October 1984.


Why do we need to increase sampling frequency at the transmitter?

Three reasons to increasing the sampling rate further are

1) To relax the requirements of the post D/A conversion filtering for image rejection.

2) Increase signal SNR by spreading quantization noise for a fixed number of DAC bits across a wider frequency range.

3) Minimize passband droop in the D/A reconstruction.

Reason 1 is the most dominant one in my opinion. The complexity and cost of a low pass filter (either digital or analog, but in this case it is the analog filter) is driven by the width of the transition band. And the reason for filtering the images is usually spectral containment as limited by regulatory requirements (and given the value of spectrum this is a HUGE reason). Reason 2 is a solid consideration if your DAC is limited by the number of bits but has wider analog BW and sufficient Spurious-Free-Dynamic Range (SFDR) to allow for overssampling benefits. Reason 3 is of least impact but I mention it as it is indeed a benefit. To address reason 3 we can often optionally simply incorporate an inverse Sinc filter digitally to compensate for the passband droop from the D/A (and this can be quite simple- see how to make CIC compensation filter).

Note interpolating DAC architectures do this for you- increase the sampling rate prior to conversion to simplify the analog anti-imaging filter requirements and give you more bits effectively since the same total quantization noise power is spread across a wider frequency range, some of which is then filtered out by the anti-imaging filter.

These graphic by Walt Kester of Analog Devices in their App Note MT-017 (https://www.analog.com/media/cn/training-seminars/tutorials/MT-017.pdf) shows this very well (Thanks Walt!):

To add, there is no benefit to the performance (sideband rejection) of the pulse shaping filter by increasing the sampling rate. The pulse shaping filters rejection is completely driven by how many symbols it spans. Two samples per symbol are all that is needed. Having more samples per symbol in pulse shaping does however simplify subsequent digital filtering if a post shaping interpolator is used- but this is exactly the same benefit described above for simplifying the analog anti-imaging filter, for the same reason.


I would say you can classify using the following main properties:

Blind Deconvolution.
Non Blind Deconvolution.

Then I'd follow:

Linear Model.

Time / Spatial Invariant Model.
Time / Spatial Variant Model.

Non Linear Model.

Time / Spatial Invariant Model.
Time / Spatial Variant Model.

I think those are the main properties of a Deconvolution Problem.
The hardest would be Blind Deconvolution with Non Linear Model of Variant Operator.


Defining Model in Kalman for Getting Position Using

In general, for position-related Kalman filters, you want your state vector to contain $x$, and $\dot{x}$ (location and velocity) components. See, for example, the Wikipedia page or this question and answer here.

If you have a measure of velocity, then it can certainly also be an input to the Kalman filter.


Noise color isn't a rigorous model for traffic noise

https://hal.archives-ouvertes.fr/file/index/docid/506615/filename/Can_2010_Applied_acoustic_71_8_p.pdf


aliasing to a 50 Hz sine wave sampled at 5000 Hz sampling rate?

You do have zeroth frequency in your signal indeed. If you don't, your signal would be centered around zero volts. Zeroth frequency is a DC component.

If you want to measure the frequency of an AC component, center the signal around zero volts (subtract the mean). By the way, FFT is not a good instrument to measure frequency of an arbitrary signal because it uses a fixed set of harmonics, which may or may not coincide with your signal (and FFT will give you a maximum on a nearest one).


Why is the median of a PDF not always equal to the mean?

If you consider the pdf as defining the distribution of mass in a thin rod of infinite length and variable density (and so the "weight per inch" of the rod is different depending on which part of the rod we are talking about), then the mean is the point about which the rod is perfectly balanced (that is, if we set up a fulcrum and place the rod on it, it will be horizontal and balanced perfectly if the fulcrum is at the mean, and will tilt downwards on one side or the other if the fulcrum is anywhere else.

The median, on the other hand, is the point for which we can say that half the mass is to the left of the median and half the mass is to the right of the median. If that is what floats your boat as the one true and the only meaning of "center of mass" and you feel that everyone else who believes otherwise is in a state of mortal sin, so be it. Just be aware that when you refer to the median as the center of mass, you are likely to be misunderstood by those lesser mortals who have not attained your state of enlightenment and who still believe that center of mass is another name for what they call the mean of the pdf or random variable.

With respect to the comments on the main question as to an asymmetric pdf means that the mean and median will be different, this answer over on stats.SE (where this question more properly belongs) gives an example of an asymmetric pdf for which the mean, median and mode all have the same value.


Your second try yielded a normalized sinc function. That function has a zero-crossing whenever $\pi(n-1)$ ia a nonzero multiple of $\pi$, which is for every $n\ne 1$ (assuming n is an integer) and is equal to 1 when $\pi(n-1)=0$.
Therefore, both solutions are the same.


It all works out because $[\mathbf{R}^{-1}]_{1, i}$ is just a number. It is the number in row $1$ column $i$ of the matrix $\mathbf{R}^{-1}$. So the product $[\mathbf{R}^{-1}]_{1, i}\mathbf{s}_i$ returns a vector of length $N \times 1$. See the line after Equation 10.


What should the bandwidth of the BPSK OFDM signal be after resampling?
  Is it 50 Hz (according to theoretical graph shown), 40 Hz (simulated
  before resampling) or 80 Hz?

What you call resampling is actually modelling D/A conversion with a discrete simulation system. Resampling consists of a repetition of every sample plus low pass filtering. Your D/A converter model makes sense but it should not be confused with resampling. The D/A converter can be described by convoluting the discrete time signal with some impulse response $g(t)$. In frequency domain this corresponds to a multiplication of $X(f)$ with $G(f)$ (where these are the Fourier transforms of $x(t), g(t)$). You chose $g(t)$ to be rectangular, so $G(f)$ is a sinc function. As $x$ is a discrete-time signal its frequency spectrum is infinite and so is $G(f)$. The first spectral nulls of $G(f)$ are at $-f_\mathrm{s}$ and $f_\mathrm{s}$. So from your plot it can be seen that a sampling frequency of 40 Hz has been assumed. The useful OFDM signal is in $[-f_\mathrm{s}/2,f_\mathrm{s}/2]$.

So to answer your question: after D/A conversion the spectrum of the OFDM signal is non-zero for arbitrarily high frequencies. To find the 3-dB bandwidth $B$ you have to find the frequency $f_0$ for which $G(f_0)=0.5$, then $B=2f_0$.

  Should a filter be used to achieve the correct bandwidth?

Yes, after D/A conversion, a lowpass filter with cut-off frequency $f_\mathrm{s}/2$ should be used to attenuate the spectral aliases that are caused by the D/A conversion. This so-called reconstruction filter is generally required after D/A conversion and is not OFDM specific. In cases where signal bandwidth doesn't play a role it can be omitted.

  How is OFDM typically implemented to ensure the correct bandwidth; should some subcarriers be set to zero? Should another DAC approach be used?

Typically a lowpass filter with cut-off frequency $f_\mathrm{c} = f_\mathrm{s}/2$ is used. Additionally, some subcarriers at high frequencies are set to zero to enable filtering at $f_\mathrm{c}$ without distorting the useful signal. This is more important at the receiver side where the bandwidth must not exceed $f_\mathrm{s}/2$ in order to avoid aliasing when sampling the received signal. An anti-aliasing lowpass filter is usually applied at the receiver to achieve that.


What is stochastic differential equation and its need?

One way to approach your problem is to use an Extended Kalman Filter with a small step size.  $x,y,$ and $z$ are your state variables.  Runge Kutta works better than Euler.  The ODE solvers in Matlab also work, but need to be set up for a single small time step.  It looks like at each time step, you propagate the differential equations forward in time, and the perturb the new states with noise, and then propagate those, repeat, ...repeat....

Google ( Lorenz Equation Matlab ) for code examples.  You have a Noisy Lorenz system.

$$ \left| \begin{array}{c} x[k+1] \\
y[k+1] \\ z[k +1] \end{array} \right|= \text{ODE STEP} \left(\left| \begin{array}{c} x[k] + \text{noise}_x \\
y[k] + \text{noise}_y \\ z[k]+ \text{noise}_z \end{array} \right|
\right)
$$

Noisey Van Der Pol Example,

This is an old homework problem with some unnecessary feature like measurement noise, but you can modify it to your needs.

There is a way for the Matlab ODE solver to only iterate from 0 to h instead of calculating 41 sub intervals.  but I don't recall how after a few years

You can use a 4th order Runge Kukta
https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

there is also
https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_method_(SDE)

You should note that, given your problem statement, you don't really have a pure SDE, The Lorenz is a chaotic system, with a noise perturbation.

function dydt = vanderpoldemo(t,y,Mu)
%VANDERPOLDEMO Defines the van der Pol equation for ODEDEMO.

% Copyright 1984-2002 The MathWorks, Inc.
% $Revision: 1.2 $  $Date: 2002/06/17 13:20:38 $

dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)];

% HW 8  part c
clear all
h=.01;
tspan = [0:h:60];
y0 = [1; 0];
epsilon = [2 1 .2] ;
q=5;
r=1;

y1=zeros(2,length(tspan));
F1=zeros(2,2,length(tspan));
y1(:,1)=y0;
ode = @(t,y) vanderpoldemo(t,y,epsilon(1));
for i=2: length(tspan)
     [t,temp] = ode45(ode, [0 h], y1(:,i-1));
    y1(:,i) = temp(end,:)'+sqrt(q)*[0 ; randn]*h;
end
y2=zeros(2,length(tspan));
y2(:,1)=y0;
ode = @(t,y) vanderpoldemo(t,y,epsilon(2));
for i=2: length(tspan)
     [t,temp] = ode45(ode, [0 h], y2(:,i-1));
    y2(:,i) = temp(end,:)'+sqrt(q)*[0 ; randn]*h;
end
y3=zeros(2,length(tspan));
y3(:,1)=y0;
ode = @(t,y) vanderpoldemo(t,y,epsilon(3));
for i=2: length(tspan)
     [t,temp] = ode45(ode, [0 h], y3(:,i-1));
    y3(:,i) = temp(end,:)'+sqrt(q)*[0 ; randn]*h;
end

figure(1)
plot(tspan,y1(1,:))
[plot_limits]=ginput(2) %select 3  orbits
startp=round(plot_limits(1,1)/h);
endp=round(plot_limits(2,1)/h);
y1=y1(:,startp:endp);
t1=tspan(startp:endp);
h1=plot(y1(1,:),y1(2,:));
set(get(h1,'Parent'),'XLim',[-4 4],'YLim',[-4,4]);
xlabel('y(1)')
ylabel('y(2)')
title([' Noisy van der Pol Equation, \epsilon = ',num2str(epsilon(1))])
%
%Plot of the second solution
 figure(2)
 plot(tspan,y2(1,:))
[plot_limits]=ginput(2) %select 3  orbits
startp=round(plot_limits(1,1)/h);
endp=round(plot_limits(2,1)/h);
 y2=y2(:,startp:endp);
t2=tspan(startp:endp);
 h2=plot(y2(1,:),y2(2,:));
set(get(h2,'Parent'),'XLim',[-4 4],'YLim',[-4,4]);
 xlabel('y(1)')
 ylabel('y(2)')
 title(['Noisy van der Pol Equation, \epsilon = ',num2str(epsilon(2))])
%
%  Plot of the third solution
figure(3)
 plot(tspan,y3(1,:))
[plot_limits]=ginput(2) %select 3  orbits
 startp=round(plot_limits(1,1)/h);
 endp=round(plot_limits(2,1)/h);
y3=y3(:,startp:endp);
 t3=tspan(startp:endp);
h3=plot(y3(1,:),y3(2,:));
set(get(h3,'Parent'),'XLim',[-4 4],'YLim',[-4,4]);
 xlabel('y(1)')
 ylabel('y(2)')
 title(['van der Pol Equation, \epsilon = ',num2str(epsilon(3))])
 z1=y1(1,:)+sqrt(r)*randn(size(y1(1,:)));  % add measurement noise
z2=y2(1,:)+sqrt(r)*randn(size(y2(1,:)));
z3=y3(1,:)+sqrt(r)*randn(size(y3(1,:)));


You may have to read on SS (spectrum sensing) techniques. There is a famous method called "Energy detection".

Else

If the noise is WGN, then its PSD is constant. So at a particular frequency of the spectrum, if the PSD is around that of noise, then it's noise, else it's signal.


(TLDR: if you want to dive in I might recommend Will Pirkle's books Designing Audio Effect Plug-Ins in C++ and Designing Software Synthesizer Plug-Ins in C++ but you might be completely overwhelmed depending on your background)

First off, if you want to create new sounds, I love Max/MSP and that will get you creating new sounds quickly but it will not teach you DSP.

Before a mod closes this as off-topic or too open-ended, I will say this:

There is no recipe, no set list of steps. Period.

Or as Umberto Eco wrote in his famous thriller "Foucault's Pendulum"....there is no map.

Take that for what you will. I'm sure everybody's story who works with audio is different. I grew up messing around with 4-track tape recorders, and since I grew up with computers, started playing with digital audio workstations like Audacity and ProTools (Ardour during my hardcore linux purist phase), making demos and playing in bands. During this time I also tried my hand at BASIC and later scripting languages. I got interested in theoretical physics in high school, and majored in physics in college (where I made a lot of demo recordings with bands and learned how to solve problems, but not a lot of domain specific knowledge (which I think is a good thing)).

I didn't want to be a particle physicist and had met a professor who worked in architectural acoustics. I found a grad program that was able to pay me to think about sound and space (mostly psychoacoustics) for 4 years and earn my MS and PhD.

At the end of my PhD I met (through a random LinkedIn message I initiated) the CTO of a startup who was willing to transfer my skills in psychoacoustics to programming microcontrollers to do effects and equalization. My job got me started down the path of really doing synthesis and processing. I had (almost) no idea what I was doing, and that's when I really learned. I'm just starting on this journey of being able to digitally manipulate sound, and I find that there's just so much out there and it is totally fascinating.

Every day I realize I know so little, but take in what I can each day (and try not to get frustrated that I don't move faster). Understanding is a slow process that cannot be expedited.

That being said, if I had to do it all over again knowing what I know now (and how could I have really?), I would have taken as much math and programming as I could have in high school (which wasn't offered in Catholic school) and majored in electrical engineering (which wasn't offered at the liberal arts college I attended). EE is really where all this stuff stems from. I'm biased but I would say it doesn't hurt to study something about hearing science and psychoacoustics too. Music appreciation can fuel your enthusiasm.


Discrete Time Fourier Transform to Frequency Domain

Most introductory texts on DSP have solved examples of this kind, so you should probably get some good text and work with it. But I'll give you a few hints anyway.

You know that

$$X(e^{j\Omega})=\sum_{n=-\infty}^{\infty}x[n]e^{-jn\Omega}\tag{1}$$

So for (a) you get

$$X(e^{j\Omega})=\sum_{n=-\infty}^{\infty}\left(\frac34\right)^n u[n-4] e^{-jn\Omega}\tag{2}$$

Since $u[n-4]$ equals zero for $n<4$, (2) simplifies to

$$X(e^{j\Omega})=\sum_{n=4}^{\infty}\left(\frac34\right)^n e^{-jn\Omega}=
\sum_{n=0}^{\infty}\left(\frac34\right)^{n+4}e^{-j(n+4)\Omega}=\\=
\left(\frac34\right)^4e^{-j4\Omega}\sum_{n=0}^{\infty}\left(\frac34\right)^{n}e^{-jn\Omega}\tag{3}$$

where you can evaluate the final sum using the formula for a geometric series with $r=\frac34e^{-j\Omega}$.

For (b) you have

$$X(e^{j\Omega})=\sum_{n=-\infty}^{\infty}a^{|n|}e^{-jn\Omega}=\sum_{n=0}^{\infty}a^ne^{-jn\Omega}+\sum_{n=-\infty}^{-1}a^{-n}e^{-jn\Omega}=\\=
\sum_{n=0}^{\infty}a^ne^{-jn\Omega}+\sum_{n=1}^{\infty}a^{n}e^{jn\Omega}=
\sum_{n=0}^{\infty}a^ne^{-jn\Omega}+ae^{j\Omega}\sum_{n=0}^{\infty}a^{n}e^{jn\Omega}$$

Now you have both sums in the desired form to apply the formula:

$$X(e^{j\Omega})=\frac{1}{1-ae^{-j\Omega}}+\frac{ae^{j\Omega}}{1-ae^{j\Omega}}$$

You can combine the two terms, which should finally result in

$$X(e^{j\Omega})=\frac{1-a^2}{1-2a\cos\Omega+a^2}$$

Please remember that this is basic stuff that you must learn to do yourself. The only way to learn it is by solving the problems which you find at the end of each chapter of any introductory DSP text.


Approximating pixel location in scale space

First, you need to ensure that the pixel has minimal value in its 8-neighborhood. SIFT, SURF and other keypoint detectors filter these pixels in step called "non-maximal suppression". This is basically a necessary condition for second-order approximation we will use to determine sub-pixel location.

If you imagine the scale space slice $I$ as a 3D surface, the pixel lays in a valley sorrounded by pixels with higer value. The shape can be well approximated by a parabola surface and the corrected pixel position is at its minimum.

The pixel value $I(p)$ at location $p=(x,y)^{\textbf{T}}$ can be approximated by second order Taylor expansion:

$$I(p+\delta)\approx I(p) + \delta^{\textbf{T}}g+\frac{1}{2}\delta^{\textbf{T}}H\delta$$

where $g$ and $H$ are gradient (2-vector of first order derivatives) and Hessian ($2\times 2$ matrix of second-order derivatives) at $(x,y)$.

The $\delta$ is a sub-pixel correction factor we want to compute.

Since the parabola has minimum at point where its derivatives are equal to zero, we obtain derivative of right-hand side of the above Taylor expansion (with respect to $\delta$) and place it equal to zero:

$$g+H\delta=0$$

From this we can compute the correction vector:

$$\delta=-H^{-1}g$$

Finally, the desired sub-pixel location is given by $p + \delta$.

The gradient and Hessian are obtained using finite differencing:

$$g_{1}=\left( I(x+1,y)-I(x-1,y) \right)\cdot 1/2$$
$$g_{2}=\left( I(x,y+1)-I(x,y-1) \right)\cdot 1/2$$
$$H_{1,1}=I(x+1,y)+I(x-1,y)-2I(x,y)$$
$$H_{1,2}=H_{2,1}=\\ \left(I(x+1,y+1)+I(x-1,y-1)+I(x+1,y-1)+I(x-1,y+1)\right)\cdot 1/4$$

So the algorithm is:

Check if the pixel is a local minimum (sorrounded by pixels with higher value).
Compute $g$ and $H$ using finite differencing formulas above.
Compute $\delta$.

You can repeat steps 2. and 3. for extra accuracy, but just one step is fine for our needs.

The value at sub-pixel location $I(p+\delta)$ can be obtained by well known bilinear or bicubic interpolation.

We did the interpolation with respect to $x,y$, but you can perform this interpolation with respect to scale as well. This will lead to 3-vector gradient and $3\times 3$ Hessian, but the Taylor expansion formula holds.

You can take a look on OpenSURF implementation in C++ or C# which has this implemented.

The above derivation is actually a Newton's Method used in multivariate optimization.