Develop a simple image classification system by using Laplacian of Gaussian (LoG)

The two simplest ideas are probably to either pad the smaller images with a constant pixel value to a common shape, or crop the bigger images to a common shape.

In the first case, you will however likely end up with a network biased towards padded images/samples with constant pixel value areas. In the second case, you loose possibly prominent information with respect to your classification task.

You could also use the cropping approach with several crops per image in order to avoid loss of information in your data set, leading to an increased proportion of the classes associated with bigger images in your training set. This also biases your neural network (this could be more or loss controlled thanks to appropriate weights attached to the classes' samples during training, if you can predict the imbalance before training, cf. this towardsdatascience article to read about training with imbalanced classes, with Keras code).

Those questions were raised and discussed in this SO post.


C++ based FIR filter design using GNU Octave remez function

There are some known bugs (such as bug #38134 and bug #48078 on savannah) with the current implementation of remez found in Octave's signal package (see remez.cc on sourceforge), not to mention that the source to link to seems to lag a little from that latest version (you are thus lacking some critical fixes that have been included to address various stability and convergence issues).

Since you appear to be reasonably comfortable with modifying sources, you may want to try out the latest remez.cc source with the submitted patch attached to bug #38134, or alternatively the newly submitted firpm replacement to remez.


Is there a way to obtain the original signal (stationary process) from its combination through filtering (matlab) and crosscorrelation?

I think that (at least in general) (and if I understand your question correctly) this is not possible. As an example, assume that $w_1(t)$ contains a harmonic component of period $t_0$. The signal $x(t)$ will not contain any trace of this component. As a result no amount of filtering $x(t)$ will be able to recover the original signal.


To get a detailed answer along the lines of what you propose, we need to be careful about the normalization used in discrete Fourier transform (DFT) and inverse discrete Fourier transform (IDFT):
$$\text{DFT: }\quad X[k] = \sum_{n=0}^{N-1} x[n] e^{-j 2\pi kn/N}\tag{1}$$
$$\text{IDFT: }\quad x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j 2 \pi k n / N}\tag{2}$$
Those have a normalization that is directly compatible with fft and ifft from MATLAB, Octave, NumPy and SciPy. The indexes $k$ and $n$ run from $0$ to $N-1$. (MATLAB and Octave have a different indexing convention, $1$ to $N$.) Then:
$$\begin{array}{c}y[n] = x[n] - x[\operatorname{mod}(n-1, N)]\\
\begin{align}\\
Y &= \operatorname{DFT}\big(y\big)\\
&= \operatorname{DFT}(x*[1, -1, 0, 0, \ldots])\\
&= \operatorname{DFT}(x)\times\operatorname{DFT}([1, -1, 0, 0, \ldots])\\
&= X\times\operatorname{DFT}([1, -1, 0, 0, \ldots])\end{align}\end{array}\tag{3}$$
$$\Rightarrow Y[k]= X[k]\left(1 - e^{-j2\pi k/N}\right),\tag{4}$$
where $\operatorname{mod}$ gives the unsigned remainder, for example $\operatorname{mod}(-1, N) = N - 1$, the symbol $*$ denotes length-$N$ circular convolution and $\times$ denotes multiplication, and all sequences are of length $N$.
Circular convolution in the discrete time domain is exactly equivalent to multiplication in the discrete frequency domain, when DFT and IDFT are used to transform the sequences between the domains. See for example Circular Convolution - MIT OpenCourseWare. For "linear" convolution of discrete sequences, there is no such elegant pair of equivalent operations, which makes me think no expression proposed as an answer to the question will ever be fully satisfactory.
Considering input signals that extend to the left of the output range $0 \le n < N$ of IDFT, calculating the backward difference using modulo indexing is conditionally equal to calculating it without it:
$$x[n] - x[\operatorname{mod}(n-1, N)] = x[n] - x[n-1]\quad\text{conditionally}\tag{5}$$
under the condition that you only calculate it for some of the indexes:
$$0 < n < N,\tag{6}$$
or for $0 \le n < N$ under the condition that there's a hint of periodicity in the signal:
$$x[-1] = x[N-1].\tag{7}$$
A sufficient but not necessary condition that satisfies Eq. 7 is that $x$ is $N$-periodic, which is however prohibited by the signal defined as finite in the question. An example of another condition that satisfies Eq. 7 is $x[-1] = x[N-1] = 0$.
@CedronDawg's first answer provides a formula which corrects in the frequency domain the error in Eq. 5 if neither condition is satisfied: $Y[k] = X[k]\left( 1 - e^{-j2\pi k/N} \right) -  x[-1] + x[N-1]$. As an alternative derivation, in length-$N$ time domain the correction is an impulse:
$$\begin{align}&\big[x[0] - x[−1] - \big(x[0] - x[N-1]\big),\, 0,\, 0,\, \ldots\big]\\
= &\big[x[N-1] - x[-1],\, 0,\, 0,\, \ldots\big],\end{align}\tag{8}$$
which in frequency domain is a constant (see DFT Pairs and Properties: pair row 2, property linearity):
$$x[N-1] - x[-1],\tag{9}$$
to be added to all elements of $Y$ calculated by Eq. 4.
For a general $x$, the condition of Eq. 6 for Eq. 5 enables to use a length $N+1$ DFT and IFT to calculate the backward difference, by shifting the input to the DFT one step to the right, and finally by shifting the output from IDFT one step to the left. With forward difference $x[n + 1] - x[n]$, the shift would not be necessary, and I think this matches your discarding of the 0th sample. For a circular convolution implementation of convolution by a finite sequence, using a longer transform is a common trick to avoid the circular effects in a sufficiently large part of the output of the IDFT. For then obtaining the DFT of a partial IDFT output, in particular a one shorter, I don't think there is any shortcut.
We could perhaps re-express the question as: What is the DFT of the length $N-1$ forward difference of $x$ of length $N$, given $x$ and $X_{N-1} = \operatorname{DFT}(x_{N-1})$, a length $N-1$ DFT of the partial sequence $x_{N-1} = \big[x[0], x[1], \ldots, x[N-2]\big]$? Analogously to Eq. 4 we have:
$$y_{N-1} = x_{N-1}*[-1, 0, 0, \ldots, 0, 0, 1]\tag{10.1}$$
$$\Leftrightarrow Y_{N-1} = X_{N-1}\times\operatorname{DFT}([-1, 0, 0, \ldots, 0, 0, 1])\tag{10.2}$$
$$\Rightarrow Y_{N-1}[k] = X_{N-1}[k]\big(e^{j2\pi k / (N - 1)} - 1\big),\tag{10.3}$$
where each sequence is of length $N-1$. The desired forward difference $f_{N-1}$ is:
$$f_{N-1} = \big[x[1] - x[0],\, x[2] - x[1],\, \ldots,\, x[N-1] - x[N-2]\big].\tag{11}$$
Eq. 10.1 can be expanded to:
$$y_{N-1} = \big[x[1] - x[0],\, x[2] - x[1],\, \ldots,\, x[0] - x[N-2]\big].\tag{12}$$
By comparing Eqs. 11 and 12, it can be seen that:
$$f_{N-1} = y_{N-1} + \big[\ldots,\, 0,\, 0,\, x[N-1] - x[0]\big],\tag{13}$$
where the sequence in brackets is of length $N-1$. Taking the DFT of both sides of Eq. 13 and applying Eq. 10.3 gives the answer:
$$\begin{align}F_{N-1}[k] &= Y_{N-1}[k] + (x[N-1] - x[0])e^{j2\pi k/(N-1)}\\
&=X_{N-1}[k]\big(e^{j2\pi k / (N - 1)} - 1\big) + (x[N-1] - x[0])e^{j2\pi k/(N-1)}.\end{align}\tag{14}$$
This is the length $N-1$ DFT of the length $N-1$ forward difference of $x$ of length $N$.
Alternatively, you might be interested in the derivative of the band-limited signal represented by the samples.


By using cyclic prefix and its associated demodulation techniques, LTE OFDM (nearly) eliminates the inter-symbol inteference (ISI).
The intra-block inteference, if I understand well your question, refers to intra-carrier inteference (ICI) when the orthogonality between sub-carriers is not maintained due to Doppler, delay spread too large, synchronisation error, ...

Ignoring equipment default that can be mitigated by using sub-carrier pilots, all remaining errors are expected to be fixed by powerful channel coding. Of course the less errors before entering to channel decoder, the better the final result is.

Please note that LTE sub-carrier spacing 15KHz and 7.5Khz are well-chosen so that it is smaller than the coherence bandwidth of typical power delay profiles. Then each sub-carrier is in flat-fading, no ICI if orthogonality is maintained.

Another thing to note that even with the latest modern channel codes, a very long codeword are required to achieve comfortable performance. Inteleaver helps to reduce the length.

So, to be short, after eliminating ISI, remaining errors are corrected mainly by channel code, with helps of pilots, inteleaver and parameter dimensionning.


A good strategy in the presence of white noise would be to do a moving average over the number of samples that equals the length of the pulse. The peak of this moving average will be your area with the noise optimally minimized.

So if you use MATLAB or Octave it would be as simple as this, assuming the pulse is 200 samples wide:

out = filter(ones(200,1),1,in);
area = max(out);

Or otherwise convolve your data with a sequence of 200 ones.

Using the Python scipy.signal package this is equivalently:

out = sig.lfilter(np.ones(200,1),1,in)
area = np.max(out)


Conversion from continuous to sampled signal?

You typically write a sine wave as
$$x(t) = sin(2 \cdot \pi \cdot f \cdot t)$$
where t is time and f the frequency of the sine wave.

Please note that both time and frequency are physical quantities, i.e. they have values and units. The units need to cancel each other since you can't apply the sine function to a unit

If you want to sample this at $f_s = 5000 Hz$ you simply calculate the value every a time interval $T_s = 1/f_s$. You get

$$x_d[n] = sin(2 \cdot \pi \cdot f \cdot n \cdot T_s) $$

or

$$x_d[n] = sin(2 \cdot \pi \cdot n \cdot f/f_s) $$

Whether this is "correct" or not depends the ratio of frequency to sample rate and whether this meets the criteria of the sampling theorem or not.


What type of information can be contained in I/Q (quadrature) data?

IQ data samples can contain double the bandwidth of information within the same frequency bandwidth and/or sample rate of single channel (non-IQ or strictly-"real") data samples.

When heterodyning a real RF signal to a strictly real signal at another frequency (usually a lower IF frequency or to baseband, for easier further filtering, processing, or demodulation, etc.), the upper and lower sidebands of the original signal around the frequency of the heterodyning local oscillator will be aliased together into half the bandwidth.  The upper and lower half of an FFT of that strictly real resulting signal will contain redundant (conjugate mirrored) result bins (only half of the full complex FFT result will contain "useful"/independant information of the aliased mix).

When heterodyning a real RF signal to an IQ signal of that same bandwidth and/or sample rate, the upper sideband and lower sideband will not be aliased together, and thus can represent twice the bandwidth. An FFT of the resulting IQ data can contain independant and different data in its lower and "upper" half of the FFT's complex result bins, representing both the unmixed upper and lower sidebands of the original RF signal.  Thus twice as much frequency information in the IQ heterodyne as compared to the single channel strictly-"real" heterodyne (of the same sample rate or result channel spectral frequency bandwidth).


Meaning of Real and Imaginary part of Fourier Transform of a signal

The real and imaginary parts of the Fourier transform of a signal $x(t)$ are the Fourier transforms of the signal's even and odd parts, respectively:

$$X_R(\omega)=\frac12[X(\omega)+X^*(\omega)]\Longleftrightarrow\frac12[x(t)+x^*(-t)]=x_e(t)\\
X_I(\omega)=\frac{1}{2j}[X(\omega)-X^*(\omega)]\Longleftrightarrow\frac{1}{2j}[x(t)-x^*(-t)]=-j\cdot x_o(t)
$$

where $X_R(\omega)$ and $X_I(\omega)$ are the real and imaginary parts of $X(\omega)$, and $x_e(t)$ and $x_o(t)$ are the even and odd parts of $x(t)$, respectively.


Acoustical Design considerations for Active Noise Cancellation Project for Air vent

I first wasn't sure that this was in the scope of a digital signal processing site, but then decided that it is about limits and environments for doing DSP, so it's quite on-topic.

Your paper (what paper? Please, please, get used to always correctly cite!) hit the spot:

  If the acoustical design of the system is not optimized,the digital controller may not be able to attenuate the undesired noise adequately.

I think you understand quite well what they're saying: You've got to deal with signals in the real world, i.e. with actual, physical, pressure waves. You're doing that by modeling the real world mathematically.

This can of course only be a success if your reality is "easy" enough to fit your model. For example, your ANC probably assumes that there's a position for the noise emitter, and position(s) for the cancellation signal emitter(s). What it probably doesn't involve is that with a bad acoustic design, noise might be leaking out of the back of your emitter, bounce off a distant wall, and come back with a huge, random phase shift.

This is science. In science, you try to reduce the amount of uncontrollable effects as far as possible to make as many things certain as possible.

Your design is a good one, in principle: Cardboard is pretty effective at absorbing sound. However, you weren't really too careful when constructing your system:

gaps. Close them.
Walls. Make them thicker.
backs of speakers, open for absorption and emission of ambient noise. Pad them.
Not too sure about the bottom of your "tube": Is it dampening cardboard or actually the table, which probably is much better at transporting power of a speaker lying flat atop than air is?
Overall system exposed to ambient noise. I'd get a big cardboard, put down a layer of insulation (e.g. styrofoam), put everything inside, fill things up with insulating material.
Is omni-directional really the pattern you're looking for in microphones? This all sounds like you want your noise source mic to be highly directional to only capture the noise, and not the effect your own emissions have at the noise source.

Then, make sure all the components, especially your sensors, work very well:

are your microphones (especially that right in front of the noise emitter) perhaps clipping?
if not 1., then, how clean is the signal your speaker produces? Maybe you want to try and produce a few tones at different volumes, and look at the spectra (coming from an RF background, I'd propose full correlation-based channel sounding, but I don't know how well that works).
Verify your model: unplug your cancellation speakers. Run your device, and record what they would have produced as sound to a file, as well as what your microphones heard. Simulate what the mics at the places where you want your noise to be cancelled should have heard based on these files. Do exactly the same with the  speakers plugged in. Compare simulation with actual microphone recording.


It's easy to show that the system defined by

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

is

Not Causal : it's because for for any value of $t$ the output $y(t)$ will always depend on future values of the input $x(t)$. Because while the dummy variable $\tau$ used in the integration ranges from $-\infty$ to $t$. The argument of the input function $x(t)$ will range from $0$ to $\infty$ due squaring. Hence for any finite value of $t$, th eoutput $y(t)$ will depend on value for $x(t)$ u to time infinity. Hence non-cusal.

Has Memory: as you guessed to due integration, current value of the output $y(t)$ depends on values other than current time $t$.

Linear: follows from the linearity of the integral operator.

Time-Varying :
Lets call shifted input as $x_2(t) = x(t-d)$ and compute the corresponding output $y_2(t)$ for the shifted input:

$$y_2(t) = \int_{-\infty}^{t} x_2(\tau^2) d\tau = \int_{-\infty}^{t} x(\tau^2 - d) d\tau$$

then also compute the shifted version of $y(t)$ by $d$ as $$y_d(t) = y(t-d) = \int_{-\infty}^{\tau = t-d} x(\tau^2) d\tau $$

Let $\beta = \tau + d$ and $\tau = \beta - d$;  substitude into the integral:
$$y(t-d) = \int_{-\infty}^{ t} x((\beta - d)^2) d\beta = \int_{-\infty}^{ t} x(\tau^2 - 2\tau d +d^2) d\tau  $$

clearly $y_2(t) \neq y(t-d)$ and system is shift-varying.

Unstable: by counter-example let $x(t) = u(t)$ and see that output goes unbounded to infinity while the input is bounded for all $t$. Hence unstable in the BIBO sense.


Is it possible to replace an integrator system with an equivalent differentiator?

Typically you can devise equivalent block diagrams in terms of just integrators or just differentiators.

A simple example is on page 2-15 of

https://web.stanford.edu/~boyd/ee102/systems.pdf

which should look remarkably familiar. You should be cautioned that mathematical equivalence doesn't mean that there aren't reasons where one form would be prefered over the other.  OP Amp integrators have better noise characteristics. You also have to band limit an analog differentiator to make it physically realizable.

Another issue is that they will not in general have the same complexity.

When modeling ordinary differential equations systems,
$$
\mathbf{M}(t) \frac{d \mathbf{y}(t)}{dt}= \mathbf{g}(t,\mathbf{y})
$$
where $\mathbf{M}(t)$ (mass matrix) is invertible is often more sparse for physical systems than
$$
\frac{d \mathbf{y}(t)}{dt}= \mathbf{h}(t,\mathbf{y})
$$
which is a more general way of writing the common state space representation.

One can write :
$$
\frac{d \mathbf{y}(t)}{dt}= \mathbf{M}^{-1}(t)\mathbf{g}(t,\mathbf{y})
$$
but if $\mathbf{M}(t)$ is sparse, $\mathbf{M}^{-1}(t)$ probably isn't ,which makes a big difference when the dimension of $\mathbf{y}$ is large.

The Matlab Control systems toolbox also includes state variable equations of the form

$$
\mathbf{E} \frac{d \mathbf{y}}{dt}= \mathbf{A}\mathbf{y}+ \mathbf{B}\mathbf{u}$$


Recognizing object of interest in object tracking using Lucas-Kanade method

Tracking objects can be simple to complex depending on what type of background you have, whether the background is static or moving, whether the object is clearly distinguishable or can share similar properties with background there by blending in the background, and whether we have unique object at hand or we are tracking multiple objects which may cross the trajectories.

The good news is, the solution to all of these exists. Depending on the nature of task at hand, you may choose simple to complex solution.

What exactly we need?

As you implement any tracking algorithm, you have an estimate of where the new position of the object is likely be. Now you need to verify that indeed the object is at the same position or different. The difference is the error that you feedback to your model. For this you need to identify/segment the object in the vicinity where your projected value. In other words you need to detect that object with in a given window.

If the object's appearance or key behaviour doesn't change quite drastically you can capture feature vectors from original reference and compare the same around the vicinity where your target is expected to be. This is search of feature vector within a window based on some distance criteria.

So the question boils down to - what should be the feature vector and what should be the distance matrix?

There are many possible answers. Writing down some of them (from simple to complex) with remark on when they might work.

Color based features: This include identifying mean or dominant color, color histogram over a local window etc. If the object color is unique and doesn't blend with background then you can identify object with this manner. See [Ref 1] and [Ref 2] below
Texture based features: There is a whole range of features here. The famous Gabor filter banks, then there is co-occurrence matrix, you can take gradient of images or find edges and have edge histograms over local regions, or you can identify Local binary patterns near the regions.
See [Ref 3]. Also, many approaches combines Color and Texture information for more effectiveness. See [Ref 4] and [Ref 5]
Shape based features: Shape based technique works when the object you are tracking has a particular shapes that distinguishes from other objects and from background. This includes many techniques. See [Ref 6], [Ref 7] and [Ref 8] below.

More complex features:

If the above features doesn't work out, specifically if the object is complex and the background as well -then we can use these features:

SIFT - The scale-invariant feature transform allows to compute key points which characterises the object and hence you can compare the same object from past frame against a new frame.  See [Ref 9]
SURF, FAST, ORB:  - If SIFT turns out to be blocked because of patent issue, and due to it's speed, then you can use it's alternatives which are SURF, FAST, ORB. See [Ref 10] for a good comparison.
HoG and Haar: Finally you can use HoG - which is a good generic object detection framework. [Ref 11] and Haar based classifiers when the problems are really complex. [Ref 12]

I am not prescribing any specific solution, and none of the above could be your direct solution - yet this will help you thinking through the problem at hand. You can pick whichever feature works for you irrespective of tracking methods that are used in those respective papers.

Go through the paper Good Features to Track  thoroughly and it will give you much better perspective on the same.

What about distance criteria?

For each different feature typically different distance measure should make sense. But, distance matrix is usually dependent on your feature vector itself.

[Ref 1]: Object Tracking Based on Color Information
Employing Particle Filter Algorithm Budi Sugandi, Hyoungseop Kim et. al.

[Ref 2]: Applying a New Spatial Color Histogram in Mean-Shift
Based Tracking Algorithm Dong Xu, Yimin Wang, and Jinwen An

[Ref 3]: Texture Features for Browsing and Retrieval of Image Data B.S. Manjunathi and W.Y. Ma

[Ref 4]: ROBUST OBJECT TRACKING USING JOINT COLOR-TEXTURE HISTOGRAM JIFENG NING, LEI ZHANG and DAVID ZHANG, CHENGKE WU

[Ref 5]: Multi-Object Tracking Using Color, Texture and Motion Valtteri Takala and Matti Pietikainen

[Ref 6] :  A Survey on Feature Extraction Techniques for Shape based Object Recognition Mitisha Patel, Purvi Tandel

[Ref 7]: RAPID AND ROBUST HUMAN DETECTION AND TRACKING BASED ON OMEGA-SHAPE
FEATURES Min Li, Zhaoxiang Zhang, Kaiqi Huang and Tieniu Tan

[Ref 8]: An Evaluation of Local Shape-Based Features for Pedestrian Detection Edgar Seemann, Bastian Leibe, Krystian Mikolajczyk, Bernt Schiele

[Ref 9]: Object Tracking using SIFT Features and Mean Shift Huiyu Zhou, Yuan Yuan, and Chunmei Shi

[Ref 10]: Image Matching Using SIFT, SURF, BRIEF and ORB: Performance Comparison for Distorted Images Ebrahim Karami, Siva Prasad, and Mohamed Shehata

[Ref 11]: Part-based Multiple-Person Tracking with Partial Occlusion Handling Guang Shu, Afshin Dehghan et. al.

[Ref 12]: Real-time hand tracking using a set of cooperative classifiers based on Haar-like features Andre L. C. Barczak and Farhad Dadgostar


Frequency of the wave in frequency domain

If you're working with a 2 dimensional data set, you need to do a 2D FFT on it (which can be dome using many 1D FFTs: doing all the rows and then all the columns, or if you prefer, the columns and then the rows).

The 2D FFT gives you a 2D array of complex numbers.  Horizontal waves will give peaks on the horizontal axis.  Vertical ones will give peaks on the vertical axis, and diagonal waves will give peaks away from both axes.


How to extract unknown packet from collision between said packet and another, known packet?

As JRE said, deconvolution is not what you need.  You need to subtract the known signal from the combined signal to leave the unknown signal by itself.

If the signals are baseband signals (I'm aware that yours are modulated- I will get to that later in the answer) then you need to know when the known signal started and how strong it is.  You can get an estimate of both of those things through cross-correlation.  You cross-correlate the combined signal with a clean reconstruction of the known signal.  You should get a strong peak at the location of the known signal.  From this you can determine when the signal started.  You can also estimate the signal amplitude using the peak's amplitude.

$$
y[n] = \sum\limits_m x[m]r[n+m]
$$

This is the definition of cross-correlation, where x[n] is the combined signal, r[n] is the reconstructed known signal, and y[n] is the cross-correlation.

$$
y[n] = \sum\limits_m (s_1[m] + s_2[m] + noise[n])r[n+m] \approx \sum\limits_m s_2[m]r[m]
$$

In this equation I have broken x[n] into its constituent parts: $s_1[n]$ (the unknown signal, $s_2[n]$ (the known signal), and noise.  If we assume that $s_1[n]$ does not correlate strongly with $s_2[m]$ then we can simply disregard $s_1[n]$ and the noise term with the understanding that they will both likely add some error to the result.

So, the peak's amplitude is approximately equal to $\sum\limits_m s_2[m]r[n+m]$.  Since you know $r[n]$ completely, you can use this value to determine the amplitude of $s_2[n]$.

Once you know the start time and amplitude of $s_2[n]$, you can subtract it out of $x[n]$, leaving only $s_1[n]$, the noise term, and whatever error was in your estimate of $s_2[n]$.

If the packets are modulated the same basic process is followed, except you will also have to find/determine any carrier offset in $s_2[n]$, and it's phase.  The phase you can get from the cross-correlation peak's phase.  The carrier offset you can find using the fourth power technique (for PSK signals), or by searching for the best cross-correlation result at various carrier offsets.


def doppler(x):
    """
    Parameters
    ----------
    x : array-like
        Domain of x is in (0,1]

    """
    if not np.all((x >= 0) & (x <= 1)):
        raise ValueError("Domain of doppler is x in (0,1]")
    return np.sqrt(x*(1-x))*np.sin((2.1*np.pi)/(x+.05))

def blocks(x):
    """
    Piecewise constant function with jumps at t.

    Constant scaler is not present in Donoho and Johnstone.
    """
    K = lambda x : (1 + np.sign(x))/2.
    t = np.array([[.1, .13, .15, .23, .25, .4, .44, .65, .76, .78, .81]]).T
    h = np.array([[4, -5, 3, -4, 5, -4.2, 2.1, 4.3, -3.1, 2.1, -4.2]]).T
    return 3.655606 * np.sum(h*K(x-t), axis=0)

def bumps(x):
    """
    A sum of bumps with locations t at the same places as jumps in blocks.
    The heights h and widths s vary and the individual bumps are of the
    form K(t) = 1/(1+|x|)**4
    """
    K = lambda x : (1. + np.abs(x)) ** -4.
    t = np.array([[.1, .13, .15, .23, .25, .4, .44, .65, .76, .78, .81]]).T
    h = np.array([[4, 5, 3, 4, 5, 4.2, 2.1, 4.3, 3.1, 2.1, 4.2]]).T
    w = np.array([[.005, .005, .006, .01, .01, .03, .01, .01, .005, .008, .005]]).T
    return np.sum(h*K((x-t)/w), axis=0)

def heavisine(x):
    """
    Sinusoid of period 1 with two jumps at t = .3 and .72
    """
    return 4 * np.sin(4*np.pi*x) - np.sign(x - .3) - np.sign(.72 - x)
In [5]:
x = np.linspace(0,1,2**11)
dop = doppler(x)
blk = blocks(x)
bmp = bumps(x)
hsin = heavisine(x)
In [6]:
fig, axes = plt.subplots(2, 2, figsize=(10,10))
ax1 = axes[0,0]
ax2 = axes[0,1]
ax3 = axes[1,0]
ax4 = axes[1,1]

ax1.plot(x,dop)
ax1.set_title("Doppler")

ax2.plot(x,blk)
ax2.set_title("Blocks")

ax3.plot(x,bmp)
ax3.set_title("Bumps")

ax4.set_title("HeaviSine")
ax4.plot(x,hsin)

for ax in fig.axes:
    ax.tick_params(labelbottom=False, labelleft=False, bottom=False,
                   top=False, left=False, right=False)

fig.tight_layout();


The way the signal $T(k)$ is generated is by applying an FIR (finite impulse response) filter to the signal $S(k)$. The transfer function of this FIR filter is

$$H(z)=1+0.1z^{-30}$$

If you want to compensate for such a filter, you need a filter with a transfer function which is the inverse of $H(z)$:

$$G(z)=1/H(z)=\frac{1}{1+0.1z^{-30}}$$

This is a recursive filter with an infinitely long impulse response (IIR). If you filter the signal $T(k)$ with this IIR filter, the result will be $S(k)$.


First of all, what do you mean when you say "Taken the first half of the resulting matrix (n/2+1)"?

The FFT will result in a vector, but maybe you mean taking only the first half of the FFT samples because the others are the negative frequencies and hence redundant?

Anyway, the process you describe is actually correct, and there is no need to do anything else. In your example with the resulting matrix a, you aparently took an FFT of size 8 (Or 16 if you then choped the second half of it). But what you do is simply take the average of all the first numbers of these vectors, and that would be an estimate of the Power spectral Density at f = 1*fs/NFFT.
Taking the average of all the second number will be an estimate of PSD at f = 2*fs/NFFT, etc.

In the example the PDS would be

PSD = 1/3*[(1+2+4) , (2+5+5) ,(3+8+6) , ... ];

When you are analizing in real time as you carry out the computation, if the signal is Wide sense stationary, then it should converge or hover around some value, but if it is not, it will deviate over time to adjust to the new coming statistics of the signal.

If i understood you correctly, that should be it.

Edit:  oops i forgot about the important scaling factor. The scaling factor should be introduced to preserve the parseval theorem between the PSD and its transform (The autocorrelation). If you windowed your signal first, i think you should apply the scaling factor 1/sum(W*fs/NFFT), where W is your window vector. Since you are merely dividing by the integral of the window function in the frequency domain This is what you add to the estimate of the frequency component  when you do the FFT by windowing first. Because the FFT alone already preserves the energy.

Edit 2:

I am even more confused with the 2-D thing. Let me try an ask you a few questions:

1 - What type of signal are you trying to process? I.e: What type of application is it? Is it a unidimensional signal, like Audio sampled at 44.1Khz? Or is it a multidimensional signal like an image or video?

2 - Where are you performing these calculations? Matlab, Octave, python, or some other?

3 - When you say 2D vector (maybe could be less confusing to call this a Matrix), What does the numbers inside mean? what are the columns and what are the rows?

4 - Did you check the link to mathworks? i think you will find what you are trying to do, is well explained there.

Also, just to make sure you understand the backgroun well enough, you could refer to some classical books that have the topic very well explained:

Digital Signal Processing, Principles, Algorithms And Applications -Proakis & Manolakis - 3rd Ed. Ch. 12

And the Father of all books in spectral Estimation:

Spectral Analysis of Signals - Stoica & Moses

Edit 3:

Ahh, ok, now i get what you are trying to do. So, here's the thing. I gave a quick read to the sections on FFT and PSD in the books by Stoica and Proakis, and still can't quite understand why they propose a scaling factor of simply 1/N to the FFT, while the matlab function for PSD using periodogram method scales the FFT by 2/(Fs*N), being N the length of the FFt, in your case 256.

So i don't know a lot of C++ but i would say do the following:

1 - Take first frame X, which consists of 256 samples of audio signal.

2 - apply 256 FFT to it. you obtain Xfft.

3 - To each complex value Xfft(i), for i = 0...255, do xfft(i)*conj(xfft(i))
which is equal to abs(xfft(i))^2. this you already did,so up to here nothing new.

4 - Multiply the whole vector Xfft by 2/(Fs*N).

This will give you the Power spectral density for your frame X. Now take the next frame and do the same, And when you reach a number of frames that you would like to plot, you stop and plot your vectors.

now since you are taking frames every 128 samples, your estimates at each time step will be sort of crammed up together, so a trick that youwould use is the fact that a Speech/audio signal is approximately wide sense statinary (WSS) in an interval of 20-50 msec (See Speech Enhancement - Jacob Benesty).
So what you could do, is assuming you have an Fs = 44100, that gives you a period of 50msec of around 2205 samples.

With the dimensions of the frames you took you could easily average 15 or 16 overlaped frames and get a Smoother estimate of the PSD. But the choice is yours, experiment with it and have fun!


Anuar Y, Welcome to the DSP community.

What you're talking about is called smoothing.
Let me explain, assume we have samples $ {\left\{ x \left[ n \right] \right\}}_{n = 0}^{N - 1} $ and we want to build estimator for $ x \left[ k \right] $ which we will define as $ \hat{x} \left[ k \right] $.

Now, we have 3 types of estimation:

The case $ k > N - 1 $ is called prediction.
This is usually what we employ Kalman Filter for.
The case $ k = N - 1 $ is called filtration.
This is easy as we basically filter data we have.
The case $ k < N - 1 $ is called smoothing.
This is basically what you're after.

So, there is plenty of information about the Kalman Filter in the Smoothing framework:

Derivation of Kalman Filtering and Smoothing Equations.
Kalman Filtering and Smoothing.
Fixed Leg Smoother.

If there is something specific you need, let us know and we'll assist you.


How can a signal have two maximum frequency components?

$x(t)$ must be a band pass signal. Under certain conditions on the sampling frequency and its relation to the lower and upper band edges of the signal, $x(t)$ can be sampled at a frequency that is lower than twice its maximum frequency, without introducing aliasing. Have a look at this article and at this answer to a related question.

For the numbers given in your question, if the lower band edge $\omega_l$ of $x(t)$ satisfies $5<\omega_l<7$ then sampling with $\omega_s=10$ will allow reconstruction of the signal without aliasing.


You might be surprised to hear that your code implements neither the conventional impulse invariance method (IIM) nor the corrected IIM. Instead it implements the matched Z-transform, which maps the poles and zeros of the analog prototype according to

$$z_k=e^{s_kT}\tag{1}$$

where $T=1/f_s$ is the sampling interval.

The two versions of the IIM can be derived as follows. If we take as an example the RIAA de-emphasis transfer function

$$H(s)=\frac{sT_2+1}{(sT_1+1)(sT_3+1)}\tag{2}$$

with $T_1=3.18\,\text{ms}$, $T_2=75\,\mu\text{s}$, and $T_3=318\,\mu\text{s}$, we first have to do a partial fraction expansion of $(2)$ to get

$$H(s)=\frac{A_1}{s-p_1}+\frac{A_2}{s-p_2}\tag{3}$$

with $p_1=-1/T_1$, $p_2=-1/T_3$, and

$$A_1=\frac{1}{T_1}\frac{T_2-T_1}{T_3-T_1},\qquad A_2=\frac{1}{T_3}\frac{T_3-T_2}{T_3-T_1}$$

With the constants used in $(3)$, the transfer function of the transformed system can directly be written as

$$H(z)=\frac{TA_1}{1-e^{p_1T}z^{-1}}+\frac{TA_2}{1-e^{p_2T}z^{-1}}\tag{4}$$

From $(4)$ it can be seen that the poles are the same as for the matched Z-transform, but the zeros are different. Note that for this example there is only one zero for the matched Z-transform as well as for the conventional IIM.

The corrected IIM according to Mecklenbräuker and Jackson subtracts a constant term from the transfer function $(4)$ of the conventional IIM:

$$H(z)=\frac{TA_1}{1-e^{p_1T}z^{-1}}+\frac{TA_2}{1-e^{p_2T}z^{-1}}-\frac{T}{2}(A_1+A_2)\tag{5}$$

Note that this solution has one more zero than the other two methods. The poles of all three methods are identical.

The figure below shows the approximation of the RIAA de-emphasis filter given by $(2)$ according to the three aforementioned transformations (IIM in red, IIM corrected in green, matched Z in magenta). IIM performs worst, matched Z-transform and the corrected IIM are similar but their errors have different signs. This fact could be used to come up with a filter that combines the corrected IIM and the matched Z-transform. Since their poles are identical, we only need to combine their numerator coefficients. The result of the combined filter is shown in black. Its approximation error is less than 1 dB over the whole frequency range.


Are there any template matching / optical flow algorithms for large image displacement?

I personally like this one. It is nicely designed for occlusions and large displacements.

But, recent trends in deep learning lead to better results in optical flow such as this and this.


What are the differences between these two Gabor filter functions

Depending on the location of the peak and the scale of the two axes of the Gaussian envelope the filter may have a large DC response. A popular approach to get a zero DC response is to subtract the output of a low-pass Gaussian filter, which is what the first of these two does. In the case of images, if the DC response is not removed, the filter will respond to the absolute intensity of the image.

This tutorial gives a bit more detail.


Using DTW on the MFCC's to match a spoken word to a template

Yes, this should be enough for a basic isolated word recognition system. Probably not something for a commercial product, but good enough for a university project or demo...
It would be better to ask the user to record a word and match against this, rather than attempt to match against a large database of utterances of the same word by different speakers. The larger the database the more likely you'd find an utterance matching the user's voice, but matching against a large database would be computationally expensive!
The first MFCC captures the signal energy - you might try to remove it when doing the matching (your feature vector will thus contain MFCC #2 to 13) - this will make recognition more robust to changes in amplitude.
You might have misunderstood something. Do you have a link to such a website?

A general proporty of MFCC is that as the coefficient index increases, the standard deviation of this coefficient decreases (the DCT can be thought of as an approximation of a PCA) - so the lower coefficients will have greater magnitude. You can compensate for that by estimating the standard deviation of each MFCC coefficient on a training set, once and for all, and dividing each coefficient by its standard deviation (or a power of it < 1.0) prior to the DTW step.


Is this piecewise-defined system time-invariant?

Let us suppose that the system with input $x$ is defined by:
\begin{equation}
  y[n] =
    \begin{cases}
      x[n] & \text{if} & n<0\\
      2x[n] & \text{if} & n\ge 0
    \end{cases}
\end{equation}
As answered by @TimWescott, the two-part  definition suggests a lack of time invariance. Finding a counter-example could be a good start. For instance a two-part signal whose behavior when multiplied by 2 is visible. Indeed a modification of the Heavide jump is interesting:
\begin{equation}
  x_0[n] =
    \begin{cases}
      -1& \text{if} & n<0\\
      \phantom{-}1 & \text{if} & n\ge 0
    \end{cases}
\end{equation}
along with the right-shifted version:
\begin{equation}
  x_1[n] =
    \begin{cases}
      -1 & \text{if} & n<1\\
      \phantom{-}1 & \text{if} & n\ge 1\,.
    \end{cases}
\end{equation}
You get:
\begin{equation}
  y_0[n] =
    \begin{cases}
      -1& \text{if} & n<0\\
      \phantom{-}2 & \text{if} & n\ge 0
    \end{cases}
\end{equation}
along with the output for the right-shifted version:
\begin{equation}
  y_1[n] =
    \begin{cases}
      -1 & \text{if} & n<0\\
      -2 & \text{if} & n=0\\
      \phantom{-}2 & \text{if} & n\ge 1\,.
    \end{cases}
\end{equation}
which are not time-shifted versions of each other.


Can amplitudes be used to differentiate one signal from another?

Taking FFT and studying the frequency of each activity would be helpful. Each activity should have its own unique frequency. You can use combination of both if any two activities have the same frequency or amplitude.


Inverse Fourier Transform From Plots (2019 edition)

That's actually a fairly a tricky one. Two ways you can go about it

Method 1

Magnitude is that of a ideal lowpass filter. Inverse of that is sinc function
Phase is that of an hilbert transformer (times -1). Impulse response is (roughly) $1/\pi t$
convolve the two

Method 2

Write out the equation for inverse Fourier transform
Pop in values from the graph and simplify a bit
You get either sum or integral for unity amplitude sine waves up to $\omega_{0}$

Neither one is pretty.


Is there some literature to understand finite state ISI channel?

you can also start by "signal processing for everyone" for Gilbert Strang, it's brief and easy to understand. Then you can go with Proakis, "Digital Communications". also there are very good online videos on youtube for India teachers, you are follow them.


Hf ionospheric channel model for large bandwidth

There is the Watterson channel model, and it can be extended to any bandwidth


We have different interleaving techniques, and matrix interleaving is one of them. But at the end all of them do one thing: interleaving is a technique to protect against burst errors (no matter how we do it).

To make it more clear, you should consider the reason a packet cannot be decoded (and is failed at the receiver). Each packet usually contains a number of codewords (for example of length $N$). These codewords are actually generated by encoding messages of length $K$ by an $[N,K]$ forward error correction (FEC) code. A given FEC  code has a limited correction capability (usually denoted by $t$) meaning that it cannot correct more than $t$ errors in each codeword.

When a packet experiences fading, a large segment of packet might get corrupted. Eventually, there might be some codewords untouched by the fading while some others are fully influenced (which is called burst error). In such case, the codewords that experience fading might have more than $t$ errors. Since the packet is only accepted if all codewords have less than $t$ errors, this results in a packet failure.

The objective by performing itrerleaving before FEC encoding and after FEC decoding is to distribute the errors evenly among different codewords. Hence, it becomes more likely for a packet to get accepted even if it experiences a "deep" fade. So Interleaving by itself is not used for error correction. It cannot change anything regarding the entropy either (since it only applies a permutation to the same random variables). It only helps to get a better efficiency from the FEC code when it faces burst errors.


Apart from a sign error, your result looks correct. The term with $(1-a_{1k})$ should have a positive sign. Also note that $\text{sgn}(x)$ as the derivative of $|x|$ is of course only valid for $x\neq 0$. If you take this into account, you can write the derivative in vector/matrix notation if you define $\text{sgn}(\mathbf{a})$ to be a vector with elements $\text{sgn}(a_i)$:

$$\nabla g=(\mathbf{I}-\mathbf{A}^T)\text{sgn}(\mathbf{x}-\mathbf{Ax})$$

where $\mathbf{I}$ is the $n\times n$ identity matrix.


I'd like to apply zero padding to it, for better frequency bin resolution.

First of all, let's state it one more time that zero padding does not improve frequency resolution of DFT. It'll only interpolate the existing spectrum on a finer frequency grid, but will not add any new information to it, otherwise.

In order to improve the true frequency resolution of a DFT computation, you have to take more samples and therefore increase the length of the signal.

Zero padding is applied by adding zeros to the end of the signal to make it length N, where you then compute $N$-point DFT of it afterwards.

Lets consider the following concrete example of the sum of two sinusoidals: $$x[n] = 0.73 \cos(0.25\pi n +\pi/12) + 0.87 \cos(0.28 \pi n + \pi/17)$$

And we shall represent $x[n]$ by a finite length $L$ segment and take $N$-point DFT of it to see a spectral vision of it.

In the first case, let signal length $L=16$ and let DFT-length be $N = 64$. The following is the resulting DFT magnitude plot:

Now, apply zero padding to $x[n]$ and make it length $N=1024$ samples (by adding zeros) and take 1024-point DFT over the still $L=16$ length signal, with the result:

As you can see although increasing DFT length (equivalently zero padding) increases the smoothness of the spectrum, it does not resolve the two sinusoidals embedded in $x[n]$ at the window length of $L=16$.

Now lets increase the signal sample count by setting $L=512$ (taking more samples of $x[n]$ instead of adding zeros to it) and again taking $N=1024$ point DFT of the segment, whose result is:

As can be seen, now the two sinusoidals are finely resolved while still computing the same $N=1024$ point DFT but over a longer signal observation length ($L=512$), hence, demonstrating the fact that the spectral resolution is improved by longer signal observation and not by longer DFT (by zero padding).


minimum-phase phase via Hilbert transform returned values

I computed and compared the minimum phase HRIR and the original one.
This is my final code:
def min_phase_coversion(HRIR):
    '''

    :param HRIR: the desired HRIR impulse response to convert into minimum phase
    :return: the minimum phase version of the original HRIR
    '''

    HRIR_fft = fft(HRIR,44100)

    #computing magnitude, tested with sinusoid, it works.
    xf = fftfreq(len(HRIR_fft), 1/44100)

    fft_magnitude = np.abs(HRIR_fft/len(HRIR_fft))

    phase=np.angle(HRIR_fft)

    minimum_phase = np.imag(-hilbert(np.log(np.abs(HRIR_fft))))

    plt.figure()
    plt.plot(xf, fft_magnitude) #computing and plotting magnitude
    plt.plot(xf, phase, label='original phase')
    plt.plot(xf,minimum_phase, label='minimum phase')
    plt.legend()
    plt.plot()
    plt.grid()
    plt.show()

    min_phase_HRIR = irfft(np.abs(HRIR_fft)*np.exp(1j*minimum_phase), 44100) # |H(w)|*e^(jPhi(w))

    HRIR_original = irfft(HRIR_fft, 44100)

    plt.figure()
    plt.plot(min_phase_HRIR[0:512], label="Min_phase", linewidth=0.5, marker='o', markersize=1) #ricostruisce bene
    plt.plot(HRIR_original[0:512], label="original phase", linewidth=0.5, marker='o', markersize=1) #ricostruisce bene
    plt.legend()
    plt.show()

    plt.figure()
    min_phase_fft = fft(min_phase_HRIR, 44100)
    min_phase_mangnitude =  np.abs(min_phase_fft/len(min_phase_fft))
    plt.plot(xf,fft_magnitude)
    plt.plot(xf,min_phase_mangnitude)
    plt.show()

    return min_phase_HRIR

The results I get seem reasonable looking at the plots:

and the magnitude is the same as I expected (the original one and the min phase one are perfectly overlapping):

However, I am not totally sure of the phase behaviour (as mentioned in my question) and its returned values.


Could someone please explain to me the relationship between frequency bands and channels?

The answer is mostly political: spectrum allocation is handled by the regional authorities based on a large number of economical, technical, political factors, and not the least on how much money can me made selling or leasing the rights. In the US that's done by the FCC.

This being said, most bands tend to have constant relative bandwidth, not absolute bandwidth. Bands that are centered in the GHz range are typically 100s of MHz wide. Bands that are centered in the 100kHZ range are maybe 10s of kHZ wide (if at all).

Using this reasoning you would expect 5 GHz band to have roughly twice the bandwidths than the 2.4 GHz band. Since the WIFI channels are constant absolute bandwidth (20 MHz), there will be more channels the higher the center frequency of the band is. But again, it's most regulatory determined.


Difference Between Iteratively Reweighted Least Squares (IRLS) and Sequential Quadratic Programming?

SQP is a method for solving smooth (objective and constraint functions are at least twice differentiable) constrained nonlinear optimization problems.  It solves a series of quadratic programming problems to converge to a solution to the Karush-Kuhn-Tucker conditions for the constrained optimization problem.

IRLS is a method for solving unconstrained minimization problems by solving a sequence of least squares problems which are obtained from the original problem by a rescaling at each iteration.


Applying zero-phase filtering in the frequency domain - also works in real-time?

It's not a loophole.  A causal linear phase filter is identical to a non-causal zero phase filter in series with a delay of half the filter length.  If you check the results of your filter with an impulse input, you will probably see such a delay.
So it's not really real-time if you count the delay, but is often close enough to seem real-time if the filter's impulse response and/or block size is short relative to your real-time requirements.


Tracking a Predefined Curve in a Series of Images

From my experience, I have successfully utilized Leo Grady's Random Walks method for this. The code is also available here. It works very well and can easily be made to run in real-time depending on the contour and image size. You could watch the video from my implementation here. Note that though, it might perfom differently than the original.


What's the rationale behind this pulse shaping function?

The pulse envelope is convolved with a Gaussian function to soften the transition between the symbols. If the standard pulse envelope were used without any softening, it would be audible to an AM detector (it runs at ~172 Hz).

Source: The IBOC Handbook: Understanding HD Radio (TM) Technology by David P. Maxson


The purpose of pulse shaping filters is not to overcome ISI as is implied in the OP's question. The only reason for using a pulse shaping filter is spectral efficiency, and in the process ISI can be introduced if not done properly. In order to limit bandwidth, pulse shaping must extend the time domain response for each pulse beyond a symbol boundary, but can be done such that the response passes through zero at the symbol decision location of ever other pulse (Nyquist shaping filter, raised-cosine filter etc).
A rectangular pulse in time is a Sinc function in frequency, with the main lobes of the Sinc that roll-off at only 1/f (20 dB/decade) as shown by the blue spectrum for a single-carrier rectangular pulse modulation in the plot below. The red spectrum shows the significant advantage is spectral containment by pulse shaping:

Notice this about the plot above to see intuitively why we are less concerned with using pulse shaping for OFDM (or when we may want to consider its use). The blue spectrum above is rolling off at 1/f where f is proportional to the symbol rate. If we reduce the symbol rate for a single carrier by a factor $N$, we equivalently increase the rate of frequency rolloff by that same factor. OFDM converts a single carrier modulation into $N$ sub-carriers, each representative of the spectrum above in blue, each with much reduced bandwidth individually and therefore overall the resulting waveform extends much less in frequency beyond its primary occupied bandwidth as given by $NR$ where $R$ is the symbol rate of each of the sub-carriers.
So to recap that with the plot of the blue spectrum above as an example: Here we see the result of a 1 MHz symbol rate (1 MSymbol/sec) single-carrier modulation, and notice within the range of the plot the spectrum rolls off to about -35 dB at 5 MHz offset. Consider if we instead sent the same amount of data by using 64 carriers, all reduced in rate by 1 MHz/64, and spaced in frequency by the same (with is orthogonal carrier spacing). The resulting spectrum would dominantly be in a 1 MHz span of frequency, with the carriers on the ends contributing the most to the spectrum going beyond that range, and specifically we would expect about -35 dBc at a 5 MHz/64 =  78.125 KHz offset. That, without any additional pulse shaping, wouldn't appear very different from the red plot above (which is a single carrier pulse-shaped waveform)! If we wanted to make the spectrum even tighter (with the diminishing return of benefit), we could certainly apply pulse shaping to reduce the adjacent leakage terms accordingly.


How to choose a phase for the deconvolution of an autocorrelation?

For a real-valued signal $C$, the autocorrelation function $M$ is a real-valued even
function and the power spectral density $m$ (Fourier transform of $M$) is a real-valued
nonnegative even function.  Now, $m(k) = c(k)c^*(k)$ where $c$ is the Fourier transform
of $C$, as you correctly assert, but given only $m$ and no other
information about $c$ (or $C$),
it is not possible to determine $c$ at all.  Note that $m(k)$ also equals $[c(k)e^{j\theta_k}][c(k)e^{j\theta_k}]^* = c(k)c^*(k) e^{j\theta_k}e^{-j\theta_k}$ for
arbitrary choice of $\theta_k$, and the inverse Fourier transform of $c(k)e^{j\theta_k}$
is also a real-valued signal if we are careful to preserve conjugate symmetry so
that we choose $\theta_{-k}$ as $-\theta_k$. In other words, there are infinitely many
quite different signals that share the same autocorrelation function, and it is not
the case that only $C$ and its time-delayed versions have autocorrelation function $M$.
As an example, all PN signals, that is, pulse trains generated from a maximal-length
linear feedback shift register (LFSR), have an inverted thumbtack autocorrelation,
and the distinct PN signals (from different LFSRs)
are even two-level signals, a property not shared
by all of the infinitely many real signals that have inverted
thumbtack autocorrelation.


Defining bound on input signal to test accumulator for BIBO stability

$B_x$ is a bound on the magnitude of the input signal. So if $x[n]=u[n]$ you have $|x[n]|\le1$ and, consequently, $B_x=1$ is a valid bound on $|x[n]|$. $B_x=0$ would mean that the input signal is zero.


Can't compute real cepstrum of real signal

It's common, when computing a cepstrum, to replace any zero's or tiny magnitudes in the 1st FFT result with some (noise) floor value to keep the scale and range of the log function "reasonable looking".

Huge negative spikes (or -inf) from the log() of tiny spectrum magnitudes don't usually provide that much added useful information to the rest of the results in cepstral-type analysis.


1-D Fourier Transform Of A 2-D Image But At An Arbitrary Orientation

You can take the 1D IFFT using a pseudo-polar grid. Fig. 2 in this paper shows an example of such grid.

Possible choices for the angle $\theta$ are going to be limited. but the transform will be free from interpolation artifacts introduced by the rotation algorithm.


Phase of the sine signal by finding FFT of it in matalb using angle

The basis function for the Fourier Transform is the complex exponential $e^{j\omega t}$
Since
$$sin(\omega t) = \frac{1}{2j}(e^{j\omega t} - e^{-j\omega t}) $$
the Fourier Transform of a sine wave a phase of -90 degrees at the positive frequency and a phase of +90 degrees at the negative frequency.
For a cosine, the phase would be 0 for both the positive and the negative frequency.


THD: estimation using periodogram or fft?

Just for answer the question, basically there is no difference at all. My problem was that I've compared a modified periodogram against a FFT. Modified periodogram does have a difference; but simple periodogram doesn't.


what is the difference between YUV format and HSV format?

So the best discussion I can point to is in this PDF on Charles Poynton's site.

The PDF goes on to much more detailed explanations of each color space, and how to convert between them. It's a bit longer than I want to type out (or cut-and-paste images of). :-)


Why is PSD estimated, and not simply computed?

Well, because in a lot of real world problems this

  If I have a deterministic signal with a fixed number of samples, shouldn't I be able to directly determine its spectral information?

is just not the case. Very often, measured signals are more of a random process. A simple and common case would be to have the desired signal and some additive noise, very often Gaussian in nature.

While you can capture the signal for a while and then deterministically calculate the PSD of those samples, what you get is the instantaneous PSD of these samples which may or may not be close to the actual PSD of the whole random process that you are sampling. It all boils down to what you want to do but in a lot of cases you want the PSD of the random process that you are sampling and also very often, just capturing some samples and calculating the instantaneous PSD just does not come close to it.

There is a lot of different methods for PSD estimation. Some may work better in some cases than other, maybe some are more computationally expensive than others. Or possibly some reach a better estimate with fewer samples as input data. In the end it boils down to select the most appropriate method for your given scenario.


Digital Audio Mixing Fundamental. What will be the data rate of mixed audio?

When you add multiple audio sources, you must first make sure that they use the same sampling rate, and then same numerical format, so that addition can be realized.
So, if they are not at the same sampling rate already, then you have to resample them to the same rate. And also make sure (by type conversions) that all of them to use the same numerical data format (same type and precision, etc.)
When you say 128 kbps alone; it looks like either an compressed data stream such as mp3, or PCM raw data at 16 kHz, 8-bits per sample. Assuming the latter is the case, then the resultant stream will again be an 16 kHz, 8-bits per sample PCM data stream, having a bit rate of the same 128 kbps.


The first component of the FFT is the DC component. The other NFFT-1 should be symmetric.
If NFFT is even, that means that

response(NFFT/2-k)=conj(response(NFFT/2+k))

for $k=1,\ldots,NFFT/2-1$ and response(NFFT/2) is arbitrary.
If NFFT is odd then

response((NFFT+1)/2-k)=conj(response((NFFT+1)/2+k+1))

for $k=0,\ldots,(NFFT-1)/2-1$.

However, MATLAB allows you to ignore all that by using the 'symmetric' argument of ifft:

response_tot = [response; zeros(NFFT-1,1)]; % assuming NFFT is odd
time_sig = ifft(response_tot,'symmetric');


What statistic is used to determine presence of a signal in noise?

This is one of the oldest signal processing problems, and a simple form is likely to be encountered in an introduction to detection theory. There are theoretical and practical approaches to solving such a problem, which may or may not overlap depending upon the specific application.

A first step toward understanding the approaches to the problem is understanding how you would measure the performance of your signal presence detector. There are two important and related metrics used to quantitatively measure how good a detector is: its probability of detection $P_d$ and its probability of false alarm $P_{fa}$.

$P_d$ is specified as the probability that your detector will indicate the presence of the signal of interest, given that the signal is actually there. Conversely, $P_{fa}$ is the probability that your detector will indicate the presence of the signal of interest, given that the signal is not there. As you might expect, then, in a perfect world, we would design a system that yields $P_d = 1$ and $P_{fa} = 0$ and call it a day. As you might also expect, it's not that easy. There is an inherent tradeoff between the two metrics; typically if you do something that improves one, you will observe some degradation in the other.

A simple example: if you are looking for the presence of a pulse against a background of noise, you might decide to set a threshold somewhere above the "typical" noise level and decide to indicate presence of the signal of interest if your detection statistic breaks above threshold. Want a really low false-alarm probability? Set the threshold high. But then, the probability of detection might decrease significantly if the elevated threshold is at or above the expected signal power level!

To visualize the $P_d$ / $P_{fa}$ relationship, the two quantities are often plotted against one another on a receiver operating characteristic curve. Here's an example from Wikipedia:

An ideal detector would have a ROC curve that hugs the top of the plot; that is, it could provide guaranteed detection for any false alarm rate. In reality, a detector will have a characteristic that looks like those plotted above; increasing the probability of detection will also increase the false alarm rate, and vice versa.

From a theoretical perspective, therefore, these types of problems boil down to selecting some balance between detection performance and false-alarm probability. How that balance is described mathematically depends upon your statistical model for the random process that the detector observes. The model will typically have two states, or hypotheses:

$$
H_0: \text{no signal is present}
$$
$$
H_1: \text{signal is present}
$$

Typically, the statistic that the detector observes would have one of two distributions, according to which hypothesis is true. The detector then applies some sort of test that is used to determine the true hypothesis and therefore whether the signal is present or not. The distributions of the detection statistic is a function of the signal model that you choose as appropriate for your application.

Common signal models are the detection of a pulse-amplitude-modulated signal against a background of additive white Gaussian noise (AWGN). While that description is somewhat specific to digital communications, many problems can be mapped to that or a similar model. Specifically, if you are looking for a constant-valued tone localized in time against a background of AWGN, and the detector observes the signal magnitude, that statistic will have a Rayleigh distribution if no tone is present and a Rician distribution if one is present.

Once a statistical model has been developed, the detector's decision rule must be specified. This can be as complicated as you desire, based on what makes sense for your application. Ideally, you would want to make a decision that is optimal in some sense, based on your knowledge of the distribution of the detection statistic under both hypotheses, the probability of each hypothesis being true, and the relative cost of being wrong about either hypothesis (which I'll talk more about in a bit). Bayesian decision theory can be used as a framework for approaching this aspect of the problem from a theoretical perspective.

In the simplest practical case, the detector might trigger a detection if the detection statistic exceeds a fixed threshold $T$. In a more complicated and practical case, the detector might have some criteria for setting an adaptive threshold $T(t)$ and trigger a detection at time $t$ if the detection statistic breaks the threshold value at that instant. In your description of your problem, you hit on one common method for setting such an adaptive threshold: calculate the neighborhood mean to estimate the "background level", then set a detection threshold some amount above that mean. This can work for some applications, and there are many other ways to arrive at such a threshold.

Given a statistical model for the detector's input and the decision rule used to map that statistic to detection conclusions, one can then calculate the detector's theoretical performance metrics. In the design phase, you would typically calculate these metrics as functions of the free design parameters that you have (for instance, the threshold $T$ above). You can then evaluate the inherent tradeoffs: "if I set $T=5$, then I get $P_d = 0.9999$, but $P_{fa} = 0.01$. That's too high of a false alarm rate, so I better increase the threshold."

Where you eventually decide to sit on the performance curve is up to you, and is an important design parameter. The right performance point to choose depends upon the relative cost of the two types of possible failures: is it worse for your detector to miss an occurrence of the signal when it happens or to register an occurrence of the signal when it hasn't happened? An example: a fictitious ballistic-missile-detector-with-automatic-strikeback-capability would be best served to have a very false alarm rate; starting a world war because of a spurious detection would be unfortunate. An example of the converse situation would be a communication receiver used for safety-of-life applications; if you want to have maximum confidence that it doesn't fail to receive any distress messages, it should have a very high probability of detection of the transmitted signal.


After upsampling with zero-insert, the higher frequency images need to be filtered out, which will serve the purpose of growing the zero insert values to the correct interpolated value in between samples. Once this is done, subsequent decimation operations can be properly done since every sample will then properly represent the signal.

For a signal that occupies BW $\in [0, \omega_1]$ with original sampling frequency at $\omega = 2\pi$ and is interpolated by $I$, after zero insert, the new sampling frequency will be at $I2\pi$ and $I-1$ images $2\omega_1$ wide will be centered at every multiple of the original sampling rate ($N2\pi$) for $N \in [1,2,\ldots I-1]$ within the new interpolated sampling rate $\omega \in [0, I2\pi]$.

In terms of the normalized angular output frequency with sampling rate $\omega = 2\pi$, the bandwidth of the original signal will be BW $\in [0, \omega_1/I]$ and the $I-1$ images will each be $2\omega_1/I$ wide and centered at multiples of $N2\pi/I$ for $N \in [1,2,\ldots I-1/I]$.

scipy.signal.firls() is a good choice for creating an interpolation filter and supports multiband filters to maximize rejection at the image frequencies. Decimation will follow a similar image folding where the $D-1$ locations that would be folded in should also be considered for filtering prior to decimation.

For further details see related post:
Choosing right cut-off frequency for a LP filter in upsampler


Creating derivative images, based on original, that can only be visually understood when stacked

...a way of taking an image (png specifically, or any other format that allows transparency) and generating N derivative images.

You can split the image into different spatial frequencies and recombine them to form the original. This is exactly the way that the Discrete Fourier Transform works.

Taking as a starting point the following image (found by Google Images and finally obtained from this article):

If you progressively filter it into different spatial frequencies, you will get intermediate images that look something like this:

Low Spatial Frequencies:

Higher Spatial Frequencies:

(...but not including the low, i.e. a bandpass filter)

Even Higher Spatial Frequencies:

(...another bandpass filter but shifted higher in the spatial frequency spectrum).

All the way to 11! Spatial Frequencies:

(...similar as above, this can be a bandpass filter or a high pass to include all remaining spatial frequencies).

Now, in your favourite image manipulation software (I am using GIMP with GMIC here), turn on transparency and put the different spatial frequency bands in different layers so that you can turn them on and off. You get something that (progressively) looks like this:

All that you need now, is some transparency "paper" and a laser printer. And also to note that if you mount the transparencies one behind the other, then you might have to adjust the scale of each image depending on the distance they will be mounted at. You can, for example, mount one on your door window and the other on the window behind the door so that when you look at the images from a specific point, they (both) make sense but the one that is further away will have to have its size slightly enlarged.

You will of course have to experiment a little bit with different spatial frequencies and different spatial frequency bands to adjust the images to the result you are trying to achieve. I only used a few bands here.

Hope this helps.

EDIT:

In the case of limited dynamic range (e.g. pixels can only be on or off), the technique remains the same but the gradient will have to be remapped through the use of halftoning.


Is digital up and down sampling linear/causal/time-invariant?

Upsampling and downsampling are both causal, since the output at time $m$ depends only on samples you've seen up time $m$.

It is also linear (apply upsampling/downsampling to a x[n]+y[n] and see that you get a upsample/downsample x[n] + upsample/downsample y[n]).

Downsampling not time invariant though - what happens if you downsample by 2 and consider inputting x[n] and x[n] delayed by 1 time unit? Now, what about upsampling?


Intuition can sometimes be misleading.  But here are some ideas that might help one move towards creating a mental picture.

An infinitely long pure sinewave in the time domain (consisting of just one frequency FT or DFT basis function) will be a single impulse in the frequency domain.

Distort the sinewave a little, but leave the waveform perfectly periodic, and the impulse will be followed by an evenly spaced harmonic series.  Usually the narrower and sharper the distortion (but keeping the waveform still perfectly periodic), the longer the harmonic series.  What might be considered a limiting case of maximum distortion, the narrowest waveform with the sharpest edge will have the longest harmonic series.  Or an infinitely long sine wave in the time domain maximally distorted into just an infinitely periodic impulse train, will produce a impulse followed by an infinitely long harmonic series, which looks a lot like another periodic impulse train.

Make it an even (cosine) function, and all the impulses in the FT will be real and thus symmetric around 0.  Add a DC offset to the distorted sine wave to complete the pulse train at 0.


I switched to using the normalized power spectrum and changed my output filtering from two-pass to one-pass. This appears to have made a massive difference:

It's difficult to tell from the picture but the response is actually much smoother and you can see that the selected bin is not jumping around in the same way. The two pass filtering I was using in particular caused a lot of ringing after sudden impulses and a one-pass filter was generally much, much better.


Help me understand how Pulse-Amplitude Modulation (PAM) works

If you define pulse amplitude modulation (PAM) as a method where a train of pulses is generated with the information coded by the amplitude of these pulses then both methods comply with that definition.

In analog PAM, the amplitudes are defined by samples of the analog source signal, whereas in digital PAM, the amplitudes are discrete and are defined by the bit pattern and some fixed mapping between bit patterns and discrete amplitude values.

In both cases, the PAM signal is given by

$$s(t)=\sum_kA_kp(t-kT)\tag{1}$$

where $A_k$ are the amplitudes representing the analog or digital source signal, $p(t)$ is the pulse, and $T$ is the symbol rate.

The signal $(1)$ is a baseband PAM signal, which can be modulated by a carrier to obtain passband PAM. Using complex symbols $A_k$ on a rectangular grid in the complex plane and two orthogonal carriers gives you quadrature amplitude modulation (QAM). Take a look at this answer for more information on passband PAM.


When I ask Matlab for a zero-padded FFT, does it (or FFTW) ever explicitly create a zero-padded input vector?

There is an FFT trick to perform a $2N$-point FFT on $2N$ real-only input samples using a single complex $N$-point FFT operation. But I've never heard of performing an $M$-point FFT on $N$ input samples without zero padding the input sequence out to length $M$. Ahmed, if you can develop a way to do that you would become famous and someone would write a song about you.


Is there any information we can extract from the residual

If your estimate is based on a linear time invariant model then you could use the following measures for the quality of your estimate: cross-correlation between the estimate and residual, and the autocorrelation of the residual. Namely the cross-correlation should be close to zero at all time, otherwise there would still be some residual dynamics that could be modeled. The autocorrelation should be close to zero at all time except at zero, where it should be close to one. If this is the case then this would mean that the residual would be white noise, which can't be predicted/modeled.


Are there analog re-configurable filters?

can we have re-configurable analog filters?

Yes. The knob you turn on your grandma's kitchen or living room radio changes the tuning of an oscillator by changing the capacitance of a component. Any stable oscillator is essentially a filter for its resonant frequency.
Also, plenty of other examples: Tunable RC filters with adjustable R; mechanically tunable metal spring filters; tunable microwave cavity resonators with tuning screws, PIN diode-based LC filters... the list goes on.
Basically, take any analog filter topology and make one of the components adjustable, you get a reconfigurable filter.
Is that useful? Well, sometimes it is, if you need to keep your passband adjustable (though that is more commonly done via octave filter banks) or need to adjust your ADC's anti-aliasing filter (though that is more commonly done using high-rate ADCs and digital filtering / sampling methodology).
Generally speaking, you avoid doing complicated things with filters in analog domain as far as possible, really. Analog filters are always a trade-off, and digital filtering hence has replaced them wherever feasible.
Can you get as much freedom as with a digital filter? Practically never. Something like a RRC is pretty trivial in digital, but will need a lot of work to even remotely approximate in analog domain, it will be dependent on temperature, and costly. The closest you'd get to an analog RRC would probably be a bucket-brigade filter, and that's almost a digital filter: discrete in time, but continuous in amplitude.


Is there a name for this smoothing formula?

1) I don't have any intuition for this smoothing formula.  This seems like some kind of windowing process that uses weights that are smaller for further away values and larger for closer values ? Is there a name for it ?

The first misconception is that this is not a smoothing formula.  It's a way of representing an already smooth continuous-time signal in discrete samples.

All it's doing is as @robert says in the comments: it's assuming that the value of the signal between $\tau_i$ and $\tau_{i+1}$ is just $Y(\tau_i)$ --- it's piecewise constant.

I believe that $p(t)$ is just a "pulse".

  2) The second equality in the smoothing relation implies that $Y(\tau_{i}) \times \sqrt  \triangle = Y_{i} $. I don't follow that ?

This is just allowing us to simplify $Y(\tau_{i}) \times \sqrt  \triangle$ and allow us to talk about $Y_i$ as shorthand notation.

edit from r b-j: it's just an alternative shorthand notation for $Y(\tau)\bigg|_{\tau = \tau_i}$ if $Y_i$ was anything else (like coefficients from some other set), i can't say that it's a ZOH with a DC gain of 0 dB.  and that "$\sqrt{\triangle}$" is chaff.

  3) It looks like the weights are going to be negative and positive depending on the value of $\tau$. Is this correct ?

The "weights" are either $0$ or $\frac{1}{\sqrt{\triangle}}$. For any value of $\tau$ outside $0 \le \tau - i\triangle \le \triangle$ the "weight" (pulse) is zero.

  4) Like I said, I can send the 3 pages ( what I wrote is a part of a longer derivation of a functional relation for least squares ) to anyone who is interested.  Dr. Kalaith makes this approximation sound trivial but it's not to me..

Please give me feedback on the above answers.


That's a (copy-paste) mistake in the book. The result in the given example is that the sum doesn't converge, yet the author concludes that the system is stable, which is obviously wrong:

$$\sum_{n=-\infty}^{\infty}|h[n]|=\sum_{n=0}^{\infty}n=\infty$$

The system is not BIBO-stable.


Converting from a phase response to a group delay response and vice-versa

The group delay of a filter is defined as minus the change in the phase response with respect to frequency.  If the phase response of a filter is $\Phi(\omega)$, the corresponding group delay $\tau_g$ is given by:

\begin{equation}
\tau_g = -\frac{d\Phi(\omega)}{d\omega}
\end{equation}

In Matlab code, the group delay of a 4th order Butterworth filter can be calculated like so:

[b a] = butter( 4, 0.25 ); % design the filter
[Hz fVec] = freqz( b, a ); % compute the frequency response and angular frequencies
phi  = unwrap(angle(Hz)); % extract phase response
taug = -diff(unwrap(angle(Hz)))/((fVec(2))); % compute group delay
figure; plot( fVec(2:end), taug );  % plot the result

The phase response must be unwrapped to avoid $\pi$ to $-\pi$ discontinuities.  Because the difference between adjacent frequency bins is the same, we need only scale the -diff(phi) by 1/fVec(2).  Finally, because the diff command produces a vector one less than its input, the group delay values are plotted against fVec(2:end) and not fVec.

Alternatively, you could use the matlab command grpdelay, though you're likely to learn more from the previous!

figure; grpdelay(b,a);


Applying Image Filtering (Circular Convolution) in Frequency Domain

In my StackExchange Signal Processing Q38542 GitHub Repository you will be able to see a code which implements 2D Circular Convolution both in Spatial and Frequency Domain.
Pay attention to the function CircularExtension2D().
This function align the axis origin between the image and the kernel before working in the Frequency Domain.
Remember that for Discrete Signals the implicit assumption on signals, In frequency Domain analysis, is being periodic (Circular).
In the discrete case one could indeed apply Circular Convolution by element wise multiplication in the Frequency Domain.
With proper padding one could apply linear convolution using circular convolution hence Linear Convolution can also be achieved using multiplication in the Frequency Domain.
See:

In depth description can be found in FFT Based 2D Cyclic Convolution.
Regarding your questions:

The filter is just an array of numbers. As long as you are after 2D Circular Convolution there is no constraints on the Filter. If it is valid for 2D Spatial Circular Convolution it is valid for Frequency Domain Circular Convolution.
Up to Floating Point Quantization errors both are mathematically equivalent (See Convolution Theorem).
If the Convolution Kernel is similar in its size to the image and both are large enough Frequency Domain Convolution becomes more efficient than Spatial Domain.


Stereo sound separation into three sources

Partially correct. Human auditory spatial perception is rather complicated. A good book on the topic would be Spatial Hearing.  How and where we localize things depends on a great many factors. For left/right localization the main effects are interaural time difference and interaural level difference. Time differences is the strongest queue. Whatever comes earlier on the left ear will be solidly localized on the left unless even if the sound is louder on the right (within reason).
If left & right signals are perfectly time aligned than your model is simplified model.


Analyze stability of a closed-loop system with Bode

I thought that, seeing the Bode plots one could tell if the closed-loop system would be stable if the 0 dB crossing occured at a lower frequency than the −180° crossing.

The usual formulation is to say that the phase margin at  the 0dB crossing is > 0. Matlab adds a vertical line to the phase diagram at the frequency where 0dB is. The length of that line is the phase margin.

The important thing to remember is that $L(s) \neq 0.1 L(s) \neq 0.01L(s)$.
By multiplying a transfer function with a number, you create a different transfer function (that describes a different system), which may or may not be stable.

The answer to the question why something multiplied by an arbitrary number, that's different from 1, does not behave as if being multiplied by 1 is because an arbitrary number different from 1 is, well, different from 1.

The phase margin is a good way to decide about stability of a system in a yes/no fashion.

But it looks like you actually want to add a variable gain $k$ to the transfer function. The question now is, for what values of $k$ is $kL(s)$ stable?

You can answer that question by choosing arbitrary values for $k$ and do the above analysis. As you experienced, this is somewhat tedious and only answering the question for discrete values of $k$.

To get the range of values for $k$, you have to approach the problem from the other direction, that is: find the gain margin from looking at the phase diagram instead of finding the phase margin from looking at the magnitude diagram.

You know that crossing the line of -180° is the critical point where the system becomes unstable at the gain of 1. Again, changing perspective now: look at the phase diagram and find the points (frequencies) where the value is -180°, now find the corresponding gains in the magnitude diagram above.

As you can see, matlab is already doing that for you. Your phase diagram intersects -180° at two points and matlab draws the gains at those frequencies into the magnitude diagram.

Changing the gain, means moving the plot in the magnitude diagram up and down. In your case, if you (choose a gain so that) move the plot down as far as the right gain margin is, the system will become unstable at that point. Moving the the plot further down (decreasing $k$ even further) will result in a an unstable system until the gain is decreased below the value of the left gain margin. The system will be stable for lower $k$ values.

I added two shifted versions of the magnitude diagram into the first body plot, one with the right gain magin (red) and one with the left gain margin (red):

If the magnitude diagram is between the red and green lines (blue area), the system is unstable. This is because the original condition for stability is false in that area (and only in that area): The phase margin at the 0dB crossing is > 0.

If you draw the root locus, you will see that there's one branch that starts in the left hand side (stable) moves into the right hand side (unstable) and later returns to the left hand side. The values for $k$ at which it crosses the imaginary axes are the two values for the gain margin.

Btw, this is why it's called "pole placement": you move the poles around in the root locus so that the system becomes stable.

tl, dr;

The bode plot has two points that are critical for stability:

the 0dB crossing in the magnitude diagram
the -180° crossing in the phase diagram

looking at either point in one diagram, and looking at the other diagram at the same frequency allows to make some assumptions about the stability of the system:

At the frequency of the 0dB crossing in the magnitude diagram, the value of the phase diagram should be > -180°. In other words: there should be a positive phase margin. The phase margin is a qualitative measure about how far away the system is from being unstable. The bigger the phase margin, the further it is away from instability.
The range of frequencies in the phase diagram with a value < -180° (In other words: a negative phase margin)  are those frequencies where the values of the magnitude diagram are not allowed to cross 0dB in order to ensure stability of the system. The boundaries of the range of frequencies (or ranges if there are many) are the -180° crossings in the phase diagram. The gain margins in the magnitude diagram at those frequencies determine the values for the boundary of the range (or boundaries of ranges) for the value of a gain $k$ for which the system will be stable.


Sketching a tri signal with horizontal shift

The peak of the tri signal is at the point where the argument of tri (that thingy inside the parentheses following tri) equals $0$. The
base of the tri signal extends from the point where the argument equals
$-\frac 12$ to the point where the argument equals $+\frac 12$ (or $-1$ to $+1$ depending on what you are calling tri).  Thus,
tri$\left(\frac{t-t_0}{T}\right)$ is an isoceles triangle whose apex
is at the point $(t_0,1)$ and whose base extends from $t_0-\frac{T}{2}$
to $t_0+\frac{T}{2}$ (or from $t_0-T$ to $t_0+T$ for the other
interpretation of tri).


phase, magnitude, audio file and information

As the comments said. "monoaural phase" is not very important and monoaural phase changes are often inaudible. "binaural" phase is very important for localisation and imaging, particularly in stereo files.

"monoaural" means the phase of a single channel or that the same phase shift or modification is applied to all channels the same way.

"binaural" means specifically the difference in the signal at the left and the right ear.


What is the real difference between DSP and AI/data science?

The purpose of both DSP and machine learning is to transform the input signal / data set into more meaningful information. This could be,

removal of noise from an audio recording
location of faces in an image
classification of objects in an image
etc

The main difference as I see it, is that in DSP the transform is designed by the engineer. The engineer will choose a set of signal processing operations that give the desired output. That choice is guided by experience and validated from the results of experiments.

In machine learning the transform is learnt. Typically, this requires a set of training signals with known outputs, over which the system will optimise its transform. The engineer still has to choose a machine learning architecture that can replicate the transform that is needed.

Often the two will be combined. And DSP is kind of machine learning anyway - the machine learning bit is done by the human and the test set is the sum of their existence.


Difficulties while understanding Convolution

The convolution sum is
\begin{equation}
y[n] = \sum_k x[k]h[n-k]
\end{equation}
If you take the term corresponding to $k=0$, the contribution is $x[0]h[n]$. Similarly the contribution of $k=1$ is $x[1]h[n-1]$. Here, $h[n]$ is delayed by $n=1$. Similarly you can see how each term of the sum contributes to overall sum. The sum is taken over values $k$ where $x[k] \neq 0$ because only those terms contribute to the sum.


How to detect the fuzzy edge in this image?

bacteria segmentation is a broad field with many approaches. here's an example. For your spesific image, I think you can do the following tricks:

d=double(imread('bac.png'));
d=d./max(d(:));

ws=7; % window size
mean_d=imfilter(d,fspecial('average',ws),'replicate'); % this is a local mean operation
d2=d-mean_d;

d3 = conv2(d2.*(d2>0), fspecial('laplacian',0),'same');
subplot(3,1,1);imagesc(d);colorbar
subplot(3,1,2);imagesc(d2);colorbar
subplot(3,1,3);imagesc(d3);colorbar

you can take it from there...


Strange noises when filtering audio signal

The context in which you use these functions is not clear, but it seems to me that your problem is "edge effects".

When you are evaluating the convolution or the biquad, you need to access samples which are outside the current buffer. Your two implementations evaluate these samples as zero. This is incorrect. For example, for the biquad, everytime you process a block of audio, you need to store the last 2 values of the input[] and output[]; and reuse them in place of input[-1], input[-2], output[-1], output[-2]. Even if the data you process comes in small chunks, you must process it as if it came in one single stream; so the state variables of your filters must not be reset to zero at the boundaries of each buffer.


Difference between step,ramp and Impulse response

Like the impulse response, the step and ramp responses describe an LTI-system. There are situations when the step response or the ramp response describe what we want to know about the system more directly than the impulse response. When the input to a system is a piece-wise constant function of time, the transitions are well described by the step response. It gives quite directly information about overshoot and settling time.


optimization of Image Reconstruction Algorithm using Genetic Algorithm in Matlab

One simple approach would be taking the mean square error (MSE) by using

fitness_1 = mean((inputimage(:) - reconstructedimage_1(:)).^2)

though, as your image size won't change, you can ommit the mean and use sum instead.

fitness_1 = sum((inputimage(:) - reconstructedimage_1(:)).^2)

In general you want to have some kind of distance measure. Look at Wikipedia (especially at Other distances) for more sophisticated approaches.


Simulation of FMCW radar in frequency domain

The function of the FMCW radar is equal to the scanning of the environment with a si-pulse. Thats why the following approach can be made:

Formulate reflection transfer function, with bandwidth $\Delta\omega$, center frequency $\omega_c$, reflection coefficient $r$, distance to target $l$, phase constant $\beta=\frac{\omega}{c_0}n_0=kn_0$ i.e.

$$R(\omega)=re^{-j2\frac{\omega}{c_0}n_0 l}= re^{-j\omega \tau n_0 }$$

Calculate the $R(\omega)$ for a wished center frequency and bandwidth in MatLab
Calculate the inverse fourier transformation of $R(\omega)$ to get $R(t)$, zero padding is recommended

You finally get the time domain response for a single target, which is a si pulse at the time $\tau=\frac{2l}{c_0}$ which can be transformed into range by the previous in line equation

Hope I could help some of those who are looking for a start.


gaussian kernel and bionomial coefficients

The well-known twentieth-century American columnist Dorothy Parker once reviewed a book
saying

  "This is not a book to be put down lightly; it should be thrown with great force."

Be that as it may, Professor Collins is suggesting that a bar plot of the binomial
coefficients on the $n$-th row of Pacal's triangle, after suitable normalization,
will look from far away like a Gaussian distribution with mean $n/2$ and standard
deviation $\sqrt{n}/2$.  He then further obfuscates this by using the fact that for
even values of $n$, we can write the standard deviation as $\sqrt{k/2}$
where $k = n/2$
while for odd $n$, it is simply $\sqrt{n}/2$, inviting the reader to arrive at
the mistaken conclusion that there is a fundamental difference between the odd
and even rows.

To get a better understanding of where the $n/2$ and $\sqrt{n}/2$
comes from, look for reading material on the DeMoivre-LaPlace approximation
which is a special case of the Central Limit Theorem.


What is the nature of the trend?

If the trend is due to DC drift (perhaps due to the signal acquisition apparatus), you may be able to remove it with a simple high-pass filter around 0.1 - 0.5 Hz. This would eliminate the need for any real-time post-processing.

DC drift gives a slowly rising/falling DC component, such as in the first figure here:


Convert image of an irregular surface to it's signal equivalent

As I can see from comments, you are taking multiple images from the same camera at diferent angle. So what you might be looking for is Structure from motion .This is similar to Stereo vision. You can find open source solutions and software tools available here. Good luck.


How do derivative masks work for finding edges in image?

Basically, the triplet 1,0,-1 (and similarly 2,0,-2) is a discrete, 3-point derivative operator. It can be applied column- (left) or row-wise (right). Let us try first with a row: a flat one, f = 0,0,0,0,0,0,0,0, a step-like s = 0,0,0,1,1,1,1,1 and one hat-like 0,0,0,1,1,0,0,0. 0 here denotes black, and 1 white.

The results will be, on the $6$ central points for each:

f: 0,0,0,0,0,0
s: 0,1,1,0,0,0
h: 0,1,1,-1,-1,0

The same reasoning applies to columns. The Sobel operators are composed of a second filter in the opposite direction: 1,2,1 that smoothes the data. So the global effect: it differentiates in one direction, smoothes in the other. THe resulting value will be high or low, positive or negative, depending on the type of contour in the image.

The (close-to-)zero values are features not detected by the derivative (hopefully non-edges). High amplitude values and their signs are hopefully the strength of the edge and the grading (positive: black-to-white, negative: white-to-black).


Scipy.signal noise with rfft compared to fft

The DFT of a real sequence is complex-valued. The array output by scipy.fftpack.rfft consists of the real part of the 0th entry, followed by the imaginary part of the 0th entry, followed by the real part of thre 1st entry, followed by the imaginary part of the 1st entry,...

Convert the output of rfft into the appropriate complex array and plot the absolute values of that array, and you will have what you seek:

yf3 = yf2[0:-2:2] + 1j*yf2[1:-1:2];
plt.plot(np.abs(yf3));
plt.show()


Which is the best method for comparing the histogram obtained after BOW?

Histogram Intersection (Pyramid Match Kernel) specifically is currently accepted as a better tool, but this doesn't generalize really. When you formulate your features or histogram computations differently, things are subject to change.


The hard one first

  What should be my sampling frequency for a given maximum 200 kHz primary signal?

Assuming your sampling is real-valued, $f_s > 2(200 + 30)\text{ kHz}$, as per the Shannon-Nyquist theorem, would sound sufficient – but the problem is that it isn't, for reasons of ambiguity, see below.

  What should be the DSP approach I should take?

First, model your resulting signal; this isn't really hard to do here. People mean different things when they say "modulation", but I'm going to go ahead and assume what seems most likely: primary signal being modulated by the secondary means that your resulting signal contains the sum frequency.

Hence, spectral estimation is really what you're after.

Furthermore, you need an estimator that is more efficient than the DFT, but doesn't have to be able to estimate a lot of spectral components – in fact, you're only interested in a single one, namely the sum of secondary and primary frequency.

I'd recommend using ESPRIT, parameterized to a single tone to detect.

Whatever spectral estimator you use, they will expose a specific variance that will typically decrease with the number of samples you take into consideration – relate your desired spectral resolution (1 Hz) to this variance, and you'll get a number of samples you'll need to get a sufficiently "exact" estimate. That number of samples will imply a sampling rate, which probably will be much higher than the aforementioned Nyquist bound.

Now, if your signal only contained the perfect primary modulated with the clean secondary signal, you'd be lost here – because if you found out your resulting signal was at let's say 150.3 kHz, you couldn't tell which part of that would be primary and which part secondary. It's mathematically impossible to tell a variation of the primary frequency from a variation of the secondary, lest your modulation doesn't just add up the frequencies. Hence, just fulfilling Nyquist for your "clean" signals won't help you – there's simply no information in these that allow you to discriminate between primary and secondary signal.

Luckily, you have noise. You don't tell us anything about that, but I'll assume the following:

It's additive to your secondary signal,
it has zero mean (which would be good from a measurement quality point of view; otherwise, we can model it as DC + zero-mean process),
it's observation's representation is symmetric in spectrum,
it's uncorrelated to both the primary and the secondary signal.

Being additive to your secondary signal implies that you'll get primary signal multiplied with noise. That means we get another, albeit "unclean" "tone"!

So we get two spectral components in our resulting signal:

A tone with frequency $f_\text{primary}+f_\text{secondary}$,
an area around the primary frequency, spectrally shaped by convolution of the spectral shape of the noise and the dirac of the primary tone

Now, 2. is a thing hard to estimate, right, because it's all "blurred out" in spectrum?

The good thing is that with sufficient oversampling and by subtracting the result of subtracting the tone estimated by ESPRIT we can find the "center" of the second component in spectrum – and that's where our "clean" primary signal was! Now, with that frequency, and our ESPRIT-estimated sum frequency, we can find the difference frequency – which is the frequency of the secondary signal.

However, your 1 Hz resolution is very unlikely to happen, for various reasons. In fact, I'm a bit too lazy to do the math now, but I don't think the time$\times$bandwidth product of your observation would allow for that. Also, the above considerations neglect any aspect of how much noise power is actually in your secondary signal, and how much additive noise you'll also see on your overall signal. My guess is that additive noise in your observation will be in the same order of magnitude as the multiplicative noise, making the estimation of a spectral center much less clear.

Also, uncorrelated multiplicative noise is usually a gross simplification – in fact, it's not how active modulators work at all.

So, really, before you try to build a system that does what you want, try to disprove the claim that you mathematically can do it at all, based on estimates of noise and signal powers, and noise bandwidths. If disproving the possibility is hard, or you get a proof that it should be possible, then go for an implementation.
If you find a proof that 1 Hz resolution is impossible, don't waste your time implementing such a system. It's natural that some phenomena aren't possible to observe with a given certainty, and that's one of the core realizations that digital communication engineers have to suffer through.


MACE correlation filter - any suggested values of of u?

The nth entry in the u vector corresponds to the value of the desired correlation peak of the nth training image. Typically you set u to be all ones if you have only images from the positive class. Entries corresponding to negative class example images are typically set to zero.