Calculate the Output of Linear Time Invariant System Given it Impulse Response

The Discrete Delta Function, $ \delta \left[ n \right] $ is the identity operator of Linear Time Invariant Systems.

Moreover, since it LTI System we can computer for each element of the filter by itself.
So the first element of the filter, $ \delta \left[ n \right] $, just outputs the signal itself.
The other element $ \delta \left[ n - 1 \right] $ just shifts the input signal.
Since the input is 1 for any $ n \geq 0 $ we subtract 1 from 1 unless it is $ n = 0 $ then we subtract zero from 1.

Hence the solution is 0, 1, 0, 0.


You usually don't buy an oscillator running at your RF – for example, there's no single "2.4 GHz oscillator" component in the wifi module of your phone. (for the low frequency of 350 MHz, you can actually build a relatively bad oscillator using a SAW filter and a transistor yourself, but don't.)

Instead, you buy some oscillator that generates a stable base oscillation, and then use a frequency synthesizer that generates a much higher frequency from that base oscillation. In the simplest case, you'd just buy a PLL chip, and a crystal oscillator running at an integer fraction of your target frequency (350 MHz).

A very much hacky way of getting a multiple of a base frequency is getting an output buffer that generates a square wave from the base oscillation, and then filter out the fundamental and all the harmonics you don't want, leaving you with just the harmonic you want. This works through the the discrete nature of the Fourier transform of the square wave. Look for the Fourier series representation of the square wave to understand where the harmonics appear!

Most RF SoCs include much more flexible designs for generation of clocks (see: fractional-N synthesis), where a much larger range of frequencies can be selectively synthesized from a reference oscillator.

The honest answer to how to build such a device is hence probably:

Look through the RF SoCs (which typically are a microcontroller, ADC/DAC and an RF frontend) of the major semiconductor producers. For example, Silabs has quite a few RF/microcontroller combos (example). That's the style of component most companies use when they need to add a sub-GHz RF interface to their devices. The example I picked can produce a whole range of frequencies, nearly continuous from 110 MHz to 965 MHz. It's not 100% an SDR, but chances are that the actual comms standard used by your remote control is indeed covered. You will, however, need to add impedance matching, antenna or connectors etc yourself. This is just a chip! Maybe there's cheap eval boards, I haven't looked into that.

Other than that: When paying for the HackRF (or other SDRs), you're paying the price for getting a device developed that works on a very large range of frequencies well, is tested, blocks out frequencies you're not interested in producing or receiving, comes with a whole firmware and driver infrastructure; when you don't need that level of sophistication, for example because you really only want to use a very limited set of frequencies and standards, things get cheaper, but you'll have to do them yourself, because they are now application-specific engineering of RF hardware.


Almost, but let's start at the beginning. If the random process $N(t)$ has a power spectrum then it is at least wide-sense stationary (WSS), i.e., its mean and its autocorrelation function do not depend on time. However, the process

$$Z(t)=N(t)\cos(2\pi f_ct)\tag{1}$$

is not stationary, and, consequently, it has no power spectrum (the mean of the random process $Z(t)$ is given by $E[N(t)\cos(2\pi f_ct)] = E[N(t)]\cos(2\pi f_c t)$ which clearly depends on time even if $E[N(t)]$ does not, a similar reasoning can be made for the autocorrelation function).

Luckily, we're usually not interested in the power spectrum of the process $(1)$, but rather in the power spectrum of the process

$$Y(t)=N(t)\cos(2\pi f_ct+\theta)\tag{2}$$

where $\theta$ is a random phase which is independent of $N(t)$ and which is uniformly distributed on $[0,2\pi)$. This random phase reflects the uncertainty of the carrier phase with respect to the process $N(t)$. So adding a random phase is not only a "trick" that works well, but it actually reflects the nature of a modulated random process much better than the process $(1)$, which assumes a known relationship between the carrier phase and the signal.

The process given by $(2)$ is WSS, and its autocorrelation function can be computed as follows:

$$\begin{align}R_{YY}(\tau)&=E[Y(t)Y(t+\tau)]\\&=E[N(t)\cos(2\pi f_ct+\theta)N(t+\tau)\cos(2\pi f_c(t+\tau)+\theta)]\\&=\frac12 E[N(t)N(t+\tau)(\cos(2\pi f_c\tau)+\cos(2\pi f_c(2t+\tau)+2\theta))]\\&=\frac12R_{NN}(\tau)\cos(2\pi f_c\tau)+\frac12 R_{NN}(\tau)\underbrace{E[\cos(2\pi f_c(2t+\tau)+2\theta)]}_{=0}\\&=\frac12R_{NN}(\tau)\cos(2\pi f_c\tau)\tag{3}\end{align}$$

where I've used the independence of $N(t)$ and $\theta$, and where $R_{NN}(\tau)=E[N(t)N(t+\tau)]$.

Finally, the power spectrum of $Y(t)$ is given by the Fourier transform of $(3)$:

$$S_{YY}(f)=\frac14\left[S_{NN}(f-f_c)+S_{NN}(f+f_c)\right]\tag{4}$$

where $S_{NN}(f)$ is the power spectrum of $N(t)$.

In sum, if we add a random phase to the carrier then the modulated process is also WSS (if the baseband process is WSS), and its power spectrum is given by $(4)$.


I think that the problem is that you are considering that $$ T = \frac{A}{1+Af} $$ is the actual transfer function of your system (an analog amplifier I guess) and this can be miss leading because that expression is just and approximation for the region of frequencys were the amplifier response is flat. If you want to make assumptions regarding the stability of the system you should work with the full Laplace transform, this is: $$ T(s) = \frac{A(s)}{1+A(s)f(s)} $$
In order to oscillate at a given frequency, the system needs to have conjugate poles on the imaginary-axis at that frequency, but in a practical case, if the system has poles on the right half-plane, as you said, the supply voltage of the amplifiers will move this poles into the j-axis

I don't have enough reputation to comment, but I think that this question that I asked a few weeks ago may be helpful.
https://electronics.stackexchange.com/questions/279365/positive-feedback-with-constant-transfer-functions


You can see this by using integration by parts:

$$\begin{align}\int_{-\infty}^{\infty}x(t)\delta'(t-T)dt&=x(t)\delta(t-T){\Big|}_{-\infty}^{\infty}-\int_{-\infty}^{\infty}x'(t)\delta(t-T)dt\\&=x(T)\delta(t-T){\Big|}_{-\infty}^{\infty}-x'(T)\\&=-x'(T)\end{align}$$

where it is assumed that $x(t)$ and $x'(t)$ are continuous at $t=T$.


Formula for calculating IIR filter coefficients for a high order low pass?

First, the equations you solve for the coefficients depend on what type of filter you're after. They also depend on how the equations are derived. For instance, a common case would be a Butterworth lowpass, converted from an analog prototype in the s-domain via the bilinear z transform. Converting a higher order lowpass with the bilinear transformation would yield a solution with more coefficients.

However, it's common practice to form higher order filters from second order sections in series, with a first order section if the filter order is odd. The responses of the individual filters multiply to give the overall higher order response. The reason for reducing a high-order filter to first and second order filter is that high-order filters are more sensitive to the limited precision of typical floating point fixed point math. The cascaded filters are all set to the same corner frequency, and the Q of each filter is set to a different values that, together, produce the intended overall Q.

Here's an article about how to produce higher-order Butterworth filters for first and second order sections, including a calculator that produces the Q values for filter orders up to 20:

Cascading filters

Here's a second-order filter calculator that you can use to generate coefficients for the individual filters:

Biquad calculator v2

Here is more detail on the bilinear z transform:

The bilinear z transform


Compensating Loudspeaker frequency response in an audio signal

Yes, you can do this with an LMS equalizer which uses the Wiener-Hopf equation to determine the least squared solution to the filter that would compensate for your channel, using the known transmit and receive sequences. The channel is the unknown being solved, and the tx and rx sequences are known.
BOTTOM LINE:
Here is the Matlab function with error checking removed:
function coeff = equalize(tx,rx,depth,ntaps)
%Determines equalizer coefficients using the Wiener-Hopf equations
%TX = Transmitted (Desired) waveform, row vector, length must be > depth+2*ntaps
%RX = Received (Distorted) waveform, row vector, length must be >=depth
%DEPTH = Depth of solution matrix (recommend 10x ntaps but based on duration of stationarity)
%NTAPS = Number of taps for equalizer filter

%force row vectors
tx= tx(:)';
rx= rx(:)';

delay=floor(ntaps/2);
A=convmtx(rx(1:depth).',ntaps);
R=A'*A;
X=[zeros(1,delay) tx(1:depth) zeros(1,ceil(ntaps/2)-1)].';
ro=A'*X;
coeff=(inv(R)*ro);

USE:
Once the coeff for the FIR filter are determined using the function above, then the Matlab filter function can process the receive sequence:
tx_recovered = filter(coeff, 1, rx)
If you want to see the channel response of the filter use:
freqz(coeff)
If you want the solution to be the estimate of the channel instead of the compensation filter that undoes the channel response, simply swap tx and rx:
coeff = equalize(rx,tx,depth,ntaps)
DETAILS FOR THE VERY INTERESTED:
See my slides below giving a high level overview / derivation of the process, This in general form is the Normal Equation (http://mathworld.wolfram.com/NormalEquation.html) used for least squared curve fitting and other applications. I believe I was first introduced to this viewpoint in demonstrating how the Normal Equation is performing deconvolution from the book "Theory and Practice of Modem Design" by John A.C. Bingham.
In practice, I typically do a cross correlation first to determine the channel response time (delay spread) and initial time alignment, and then use an initial equalizer FIR length (# of taps) that exceeds the delay spread (not knowing if leading or trailing echos dominate I will typically start with 2x the delay spread for the FIR length). Once I see the result, the filter size can be reduced if desired based on insignificant magnitudes of the coefficients at the edges of the filter.  If the sequences are not exactly aligned, but still within the span of the filter, then the dominant tap will be offset accordingly- so not critical to align beforehand and this gives you insight into what happens if they are grossly misaligned.

Here is an interesting example of the equalizer function I used recently on a sound file from Dalen to equalize the waveforms received by the left and right channels as received by two microphones (treating left as transmit and right as receive and ignoring the actual third party transmitter for the two). The two channels are not recognizable prior to equalization, and completely aligned in amplitude, phase and characteristic after.
Here is a plot of the left and right channels prior to equalization:

Here is the same plot after equalization, right was filtered with the equalizer, and left was filtered with a simple filter just as long as the equalizer with a single unity gain tap in the center and zero elsewhere (to match the delay as the equalizer assumes nominal delay is in the center of the equalizer filter):

This is a zoom in plot of the waveforms after equalization showing how identical the two sequences have become:


What is the frequency transfer function when $x(t) = e^{-3t}u(t)$ and $y(t) = 2u(t)[e^{-t}-e^{-4t}]$

Your answer is correct; it's just in a different form than the given answer. In order to get your answer in the given form, do the following:

rewrite $Y(j\omega)$ by combining the two terms: $Y(j\omega)=\displaystyle\frac{N(j\omega)}{(j\omega+1)(j\omega+4)}$, where $N(j\omega)$ is a (very simple) polynomial. I'm sure you know how to obtain $N(j\omega)$.
write the frequency response as $\displaystyle H(j\omega)=\frac{N(j\omega)(j\omega+3)}{(j\omega+1)(j\omega+4)}$ and use partial fraction expansion to obtain the answer in the given form.

Alternatively, you could also determine $h(t)$ via inverse Fourier transform of $H(j\omega)$ as you obtained it (and noting that multiplication by $j\omega$ corresponds to differentiation in the time domain). You will see that $h(t)$ is in the form $A_1e^{-t}u(t)+A_2e^{-4t}u(t)$, which directly leads to $H(j\omega)$ in the given form.


Yes, due to the separability of the kernel :

$$e^{-j \left(\frac{2\pi}{N_1} n_1 k_1 + \frac{2\pi}{N_2} n_2 k_2 \right) } = e^{-j \frac{2\pi}{N_1} n_1 k_1} \cdot e^{-j \frac{2\pi}{N_2} n_2 k_2}$$

the 2D-DFT sum
$$X[k_1,k_2] = \sum_{n_1} \sum_{n_2} f[n_1,n_2] e^{-j \frac{2\pi}{N_1} n_1 k_1} e^{-j \frac{2\pi}{N_2} n_2 k_2} $$

can be implemented using row-column (or column-row) decompositions.

$$X[k_1,k_2] = \sum_{n_1} \left( \sum_{n_2} f[n_1,n_2] e^{-j \frac{2\pi}{N_2} n_2 k_2} \right) e^{-j \frac{2\pi}{N_1} n_1 k_1} $$

The following MATLAB/OCTAVE code shows how to apply :

N = 8;            % length of columns-rows
x = randn(N,N);   % row data

X = fft(x,N) ;    % 1D-fft along columns of x[n1,n2]
X = fft(X.',N).' ;% 1D-fft along rows of intermediate X[k1,k2]


Why can't music recognition services like Shazam or SoundHound recognize music if its speed is changed by more than 5 percent?

this is a guess from someone who does audio, but has never ever done recognition nor watermarking algorithms for someone in the industry (but i have a little idea how they may be done).

my guess is that they extract some features from the STFT of the music.  associated with each feature is an envelope that represents the quantity of that feature at the center window time of the STFT.

the collection of those envelopes can be thought of as a point in a metric space.  if you have a metric space, you can define a difference function

$$ \| \mathbf{X} - \mathbf{X}_n \| $$

between those of any pair.

then they compare the test envelope, $\mathbf{X}$, with each template, $\mathbf{X}_n$, of a set that come from the original targets (the copyrighted songs they're fishing for).  there is a norm taken of the difference between these points in a metric space, between the collection of envelopes taken from the test sound and the template.  if the norm is below a given threshold and if the norm is the minimum compared to all of the templates, then the song of that best matching template is identified as matching the input or test sound.

now if the pitch is wrong it won't match as good.  but they could make their difference algorithm insensitive to a constant offset to pitch.  it's like running the difference signal through a DC blocking filter before calculating the norm of it.  but they would have to do that, otherwise a constant pitch difference between test tune and target will always be an error.

and they could match a tempo difference, but that is much harder.  if memory and computational effort is not a problem (and they often are not for non-realtime applications running on a GHz-speed computer or server) then there could be a template for each of many speeds.  and see which one matches best.  but if these other tempo templates do not exist, and you submit a tune that does not exactly match a template of it at the original speed, then at some point in the tune or song, the pitches will not match at all.  that will make a dumb algorithm think they're different tunes.


Phase Modulation and Intermodulation Anomaly ? (Spectrum plots inside - **Updated**)

A couple of points to think about:

1) Looking at the scales of your FFT plots, it appears that the fundamental is around -25 dBW, whereas the 5th harmonic is around -125 dBW. This is a very large difference, implying the system is only weakly nonlinear. In practice you'd have problems detecting the 5th harmonic even without the phase modification just due to thermal noise. If you change your scales to have a dynamic range of say 100 dB, you would see the two plots look nearly identical.

2) The phenomenon you are modelling looks to be AM-PM conversion. Under this interpretation, it can be thought of as "adding" some phase noise to the output signal. (The tails of the phase noise spectrum would droop down and mask weaker signals.) This could be an explanation for why the result looks like a noise floor.

3) A sanity check could be-- what does the spectrum look like when the peak phase shift is halved? Or cut to 1/4th, etc? Does the noise floor drop by some commensurate amount?


How to calculate FFT of a non-uniformly sampled signal?

You are essentially referring to nonuniform signal sampling in which data samples are not acquired at exact periods but have same random jitter deviations from their exact uniform timing positions.

Now under some suitable conditions (such as the jitter length being less than a period), there are very effective algorithms which can perfectly convert the nonuniformly sampled data into their equivalent uniform samples. (i.e., given the values of nonuniform samples they give a new set of values which correspond to the perfect uniform samples).

After this nonuiform to uniform conversion, you can then apply the usual FFT, which assumes that the data to be transformed was uniformly sampled.

However, before such a DSP attempt, I would suggest you understand the sampling and read out mechanism of the Arduino Uno and MPU 6050 6DOF sensor combination. That sensor can be read in either a single data read mode or in a burst buffer data read mode. In this second mode, I think it will provide a set of buffered data which are sampled by the MPU6050 unit itself as uniform as possible, compared to a sigle data read mode mastered by the Arduino.

From your data, there seems to be a drift of a few milliseconds per sample. I assume a 20 milli second of sampling period is quite achievable with less than 1 ms of timing accuracy using the Arduino in a proper configuration. So may be your main problem is on the software configuration?

EDIT: Try the following MATLAB code as is:
NOTE: The following code MAY result in WORSE results than available before processing... The success highly depends on the application.

% NONUNIFORM SIGNAL DEFINITION
% ----------------------------
N=100;       % Converts a block of length N nonuniform samples to uniform
K=N;        %
M = N*K;    %
Tsam = 0.020; % Sampling period in seconds,

% tp represents the nonuniform sampling times
tp = [0:Tsam:(K-1)*Tsam] + Tsam*(rand(1,K)-0.5);
xnunds = sin(2*pi*1*tp);        % A simulated nonuniformly sampled data...

% GENERATE THE MATRIX AND ITS INVERSE
% -----------------------------------
A=zeros(K,K);       % initial matrix to be inverted
for p=1:K
    for q=1:K

        A(p,q)=sinc((tp(p)-tp(q))/Tsam);

    end
end
AI=inv(A);          % Inverse matrix is here

% APPLY THE CONVERSION FORMULA
% -----------------------------
xint=zeros(1,K);    % interpolate to uniform samples

for n=1:N:M

    ACS=0;

    for m=1:K       % per nonuniform sample loop

        PSI=0;

        for q=1:K   % calculation of PSI for each m

            PSI = PSI + AI(q,m) * sinc( ((n-1)*Tsam/N - tp(q))/Tsam);

        end%q

        ACS=ACS+xnunds(m)*PSI;  % acumulation of multiplicant

    end%m

    xint( (n-1)/N + 1)=ACS;    % result for each n is assigned
end%n

figure,stem(xnunds);
hold on
stem(sin(2*pi*1*[0:K-1]*Tsam),'g');
title('nonuniform samples imposed on true samples');

figure,stem(xint,'r');
hold on;
stem(sin(2*pi*1*[0:K-1]*Tsam),'g');
title('interpolated uniform samples imposed on true uniform samples)


Why we add imaginary part in inverse fourier transform since the time domain signal has only real value?

then the imaginary part must go to zero

No. That's a fallacy. There's really not much to explain here other than this intuition is wrong.
The imaginary part of a Fourier transform is only zero if the transformed signal has conjugate symmetry (which is the same as symmetry for real signals) around t=0.
Think about it. You decomposed your $e^{-ix}=\cos x + i\sin x$, and then used an integral to project your signal onto these two things. (since the integral of a product between a signal and another (conjugated) signal is the inner product of the signal vector space, and hence the projection, if you're after the mathematical formalism.)
The part of your signal that is symmetrical to t=0 projects onto the cosine, and the part that is antisymmetrical to t=0 to the sine. The sine has an imaginary unit in front.
This is really among the first things that are typically illustrated in textbooks when they introduce the Fourier transform, so this might be a great point to go back to that!


Improving Template Matching Algorithm for OpenCV

I assume that the image have same scales in terms of sizes.
With that out of the way, I would suggest try using the following:

Break up your image into parts. (the template image as well as the test image. may be 3 or 4 according to your choice.)
Then perform template matching on each part. Consider common threshold for every template part and then count it as positive only when all parts give a positive result.

(P.S. A similar solution was suggested on some other post. I don't exactly remember which post.)

Update: Partition your template and test images at exactly same points.


The $u(t-2)$ term needs to be included because the upper limit of the definite integral is $t-2$, not $\infty$:

\begin{align*}
h(t) &= e^{-(t-2)}u(t-2)\\
y(t) &= x(t) * h(t)= \int_{-\infty}^{\infty} x(t')h(t-t') dt'\\
 &= \int_{-\infty}^{\infty} x(t')e^{-(t-t'-2)}u(t-t'-2) dt'\\
u(t-t'-2)&= \begin{cases} 1 & t'<t-2 \\ 0 & \text{otherwise} \end{cases}\\
y(t) &= \int_{-\infty}^{t-2} x(t')e^{-(t-t'-2)}(1) dt' + \int_{t-2}^{\infty} x(t')e^{-(t-t'-2)} (0) dt'\\
&= \int_{-\infty}^{t-2} x(t')e^{-(t-t'-2)} dt'.
\end{align*}


Compressive Sensing with Square Measurement Matrix

It's important to clearly distinguish the terms Compressive Sensing (CS) and Sparse Signal Recovery (SSR). CS is about taking fewer measurements than what classical criteria such as Nyquist would dictate us. In the discrete setting this means taking $M<N$ measurements of an $N$-sparse vector.

SSR is about how to recover our $N$-sparse vector from your measurements.

That said, this should answer your question: If you apply CS, then you need SSR to recover your signal, however, the converse is not true - you can apply SSR even if you didn't do any subsampling. For SSR, your measurement matrix can be flat or square, it could even be tall.

Mathematically speaking, for SSR you solve a problem of the sort $$\min_x \|y - A x \| + \lambda h(x),$$ where $h(x)$ is a regularizer that encourages sparse solutions, i.e., an $\ell_0$-"norm" or its relaxation, the $\ell_1$-norm. In CS, such problems are solved since $A$ is flat and hence we must regularize to find our sparse solution out of the infinitely many possible solutions. The same problem can be used in a non-CS setup where $A$ is square or even tall, e..g, to find an approximate solution that shows a desired level of sparsity.


Jumping from DSP concepts to practice...recommend books?

The USRP with Gnu Radio software is a good and relatively inexpensive (not that cheap, but much cheaper than almost anything else out there) DSP platform.  Gnu Radio is also open source, so if you develop something good you can submit it to the project.


Find tracing/plotting stacked bands/lines in an image

Why not try just summing the columns of the image?

Well, you could try smoothing the image first, but that will have the same effect as just smoothing the output.

R code for this is below the first lot.

Here's the extra image, processed.

R Code Below

# 26593
# install.packages("png")

file <- "Q26593_original.png"

library(png)

img <- readPNG(file)

d <- dim(img)

par(mfrow=c(2,1))
plot(1:2, type='n', main="Plotting Over an Image", xlab="x", ylab="y")
lim <- par()
rasterImage(img, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])

plot(colSums(img[,,1]), type="l")

R Code #2 Below

# 26593
# install.packages("png")
# install.packages("EBImage") --- Didn't work.
# ## try http if https is not available
# source("http://bioconductor.org/biocLite.R")
# biocLite("EBImage")

file <- "Q26593_original.png"

library(png)
library(EBImage)

img <- readPNG(file)

d <- dim(img)

#xx <- (0:(d[2]-1) - d[2]/2)/d[2]
#yy <- (0:(d[1]-1) - d[1]/2)/d[1]
#sigma <- 0.2

#w <-  exp(-yy*yy/sigma) %*%  t(exp(-xx*xx/sigma))
#w1 <- exp(-yy*yy/sigma) %*% t(rep(1,d[2]))

f <- makeBrush(11, shape='disc', step=FALSE)
imgf <- filter2(img, f)
imgf <- imgf/max(imgf)

par(mfrow=c(3,1))
plot(1:2, type='n', main="Filtered Image", xlab="x", ylab="y")
lim <- par()
rasterImage(imgf, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])

#plot(1:2, type='n', main="Window Function", xlab="x", ylab="y")
#lim <- par()
#rasterImage(w1, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])

plot(colSums(imgf[,,1]), type="l", main="Column sum of filtered image")

plot(filter(colSums(img[,,1]), rep(1,11)), type="l", main="Column sum of unfiltered image, filtered")


When choosing one of these 4 types of linear phase filters there are mainly 3 things to consider:

constraints on the zeros of $H(z)$ at $z=1$ and $z=-1$
integer/non-integer group delay
phase shift (apart from the linear phase)

For type I filters (odd number of taps, even symmetry) there are no constraints on the zeros at $z=1$ and $z=-1$, the phase shift is zero (apart from the linear phase), and the group delay is an integer value.

Type II filters (even number of taps, even symmetry) always have a zero at $z=-1$ (i.e., half the sampling frequency), they have a zero phase shift, and they have a non-integer group delay.

Type III filters (odd number of taps, odd symmetry) always have zeros at $z=1$ and $z=-1$ (i.e. at $f=0$ and $f=f_s/2$), they have a 90 degrees phase shift, and an integer group delay.

Type IV filters (even number of taps, odd symmetry) always have a zero at $z=1$, a phase shift of 90 degrees, and a non-integer group delay.

This implies (among other things) the following:

Type I filters are pretty universal, but they cannot be used whenever a 90 degrees phase shift is necessary, e.g. for differentiators or Hilbert transformers.
Type II filters would normally not be used for high pass or band stop filters, due to the zero at $z=-1$, i.e. at $f=f_s/2$. Neither can they be used for applications where a 90 degrees phase shift is necessary.
Type III filters cannot be used for standard frequency selective filters because in these cases the 90 degrees phase shift is usually undesirable. For Hilbert transformers, type III filters have a relatively bad magnitude approximation at very low and very high frequencies due to the zeros at $z=1$ and $z=-1$. On the other hand, a type III Hilbert transformer can be implemented more efficiently than a type IV Hilbert transformer because in this case every other tap is zero.
Type IV filters cannot be used for standard frequency selective filters, for the same reasons as type III filters. They are well suited for differentiators and Hilbert transformers, and their magnitude approximation is usually better because, unlike type III filters, they have no zero at $z=-1$.
In some applications an integer group delay is desirable. In these cases type I or type III filters are preferred.


Amplitude distribution of Rician faded complex signal

This was a scaling problem. I did not have the same vectors for X-axis, thus creating uneven areas below PDF curves.


Cross power spectral matrix from PSDs and coherence using matrix notation

$S= diag(\underline{s}) \cdot \underline{\underline{\gamma}} \cdot diag(\underline{s})$


How can I graph the amplitude for a given frequency window over time?

Do you have access to MATLAB or Octave or similar?

In general, I would:
 - Compute the spectrogram (STFT) of the signal you're interested in.
 - Perform any required scaling to get a good amplitude estimate.
 - Apply some thresholding to treat the noise floor.
 - Devise a scheme for identifying the amplitude maxima of the harmonics you're interested in.

This final point will be much easier if you know the fundamental. For example, if you're playing a note at 440Hz (A440), and you know you only want to examine the fundamental, second, and third harmonics, you can easily define some brackets to search within.

If you don't know the fundamental, you will have to implement some method for determining what the fundamental is. However, assuming you do know the fundamental, I've provided a MATLAB example.

This example first builds an arbitrary signal with some sinusoids that decay at differing rates, then analyses the signal to identify harmonic amplitudes over time.

%% Generate a signal with harmonics.
% Build fundamental and harmonic signals.
samplePeriod = 0.001;
time = 0:samplePeriod:5; % [s]
harmonics = [1, 3, 5, 7];
fundamentalFreq = 50; % [Hz]
componentFreqs = fundamentalFreq .* harmonics; % [Hz]
componentFreqs = componentFreqs .* 2 .* pi; % [rad/s]
componentFreqs = repmat(componentFreqs', 1, numel(time));
componentSignals = sin(componentFreqs .* repmat(time, numel(harmonics), 1));

% Apply decay to the harmonics.
decay = [0, 0.1, 0.5, 1];
harmonicAmplitude = zeros(numel(harmonics), numel(time));
for iHarmonic = 1 : numel(harmonics)
    harmonicAmplitude(iHarmonic, :) = exp(-1 * decay(iHarmonic) * time);
    componentSignals(iHarmonic, :) = componentSignals(iHarmonic, :) .* harmonicAmplitude(iHarmonic, :) .* 1;
end

% Generate final signal.
signal = sum(componentSignals, 1);

%% Perform time-frequency analysis.
sampleFrequency = 1 / samplePeriod;

% Compute the spectrogram.
NFFT = 2^12;
winDuration = 0.1; % [s]
winLength = winDuration * sampleFrequency;
winFactor = mean(hamming(winLength));
[S, freq, specTime, p] = spectrogram(signal, winLength, [], NFFT, sampleFrequency);

% Compute FFT amplitude.
amplitude = abs(S) .* 2 ./ winLength ./ winFactor;

% Threshold.
floorLevel = max(max(amplitude)) * 0.1;
amplitude(amplitude < floorLevel) = 0;

% Determine the amplitude of the harmonics of interest.
fundamentalFreq = 50; % [Repeat of above]
harmonicsOfInterest = 1 : 10;
harmonicSearchWidth = 5; % [Hz]
harmonicBracketLower = harmonicsOfInterest .* fundamentalFreq - harmonicSearchWidth / 2;
harmonicBracketUpper = harmonicsOfInterest .* fundamentalFreq + harmonicSearchWidth / 2;

% Find the amplitude maxima for each harmonic at each point in time.
estimatedHarmonicAmplitude = zeros(numel(specTime), numel(harmonicsOfInterest));
for iHarmonic = 1 : numel(harmonicsOfInterest)
    activeBracket = freq >= harmonicBracketLower(iHarmonic) & freq <= harmonicBracketUpper(iHarmonic);
    bracketAmplitude = amplitude(activeBracket, :);
    estimatedHarmonicAmplitude(:, iHarmonic) = max(bracketAmplitude, [], 1)';
end

% Visualise output.
figure;
surf(specTime, freq, amplitude, 'EdgeColor', 'none');
xlabel('Time [s]');
ylabel('Frequency [Hz]');
zlabel('Amplitude');
view(2);

figure;
plot(specTime, estimatedHarmonicAmplitude);
xlabel('Time [s]');
ylabel('Amplitude');
box on;
grid on;
legendBuilder = @(a) ['Harmonic: ', num2str(a, '%1.0f')];
for iHarmonic = 1 : numel(harmonicsOfInterest)
    legendString{iHarmonic} = legendBuilder(harmonicsOfInterest(iHarmonic));
end
legend(legendString);


The best way to compare your method to others is to run the algorithms on a common data set. Quntative measures such as PSNR vary significantly across images.

It is also important to consider subjective analysis. After all, your research is aiming for improvements in appearance.


Recover signal from its power spectral density

Generally (in the absence of prior knowledge) , I don't think this could be done.

The PSD is the absolute value  squared of the DFT of the signal, and the abs function is not 1-to-1. For example abs(x) could be produced by f=x and f=-x.

However, if you have some prior knowledge about the signal then you could reconstruct it from the PSD. Without more knowledge of what you're trying to do I don't think I could help you.


The physical meaning of comparing the phases of two signals with *different frequency*?

There is not a phase in the usual sense in which there is a number $\phi_0$ such that

$$x_1(t)=x_2(t-\phi_0) \ \forall t$$

where $x_1(t)$ and $x_2(t)$ are two periodic signals of different frequencies.

However, there are some phase-like analysis one could perform, depending on the application. I'll talk about two examples that come across my mind quickly.

Suppose we have the two signals $x_1(t)=\sin(t)$ and $x_2(t)=\sin(3t-\phi_0)$. As the latter is the 3rd harmonic of the former, their sum is periodic. However, the shape of the resulting wave would vary a lot depending on $\phi_0$ (some kind of "phase difference" between the signals). In the image below, I plotted three different waves ($x_1(t)+x_2(t)$ for three different values of $\phi_0$):

$\phi_0=0$ in red.
$\phi_0=\pi$ in blue.
$\phi_0=\pi/2$ in green.

As you can see, there is a strong dependence on the value of $\phi_0$. This would be very important in music, for example, as the resulting wave would have a pretty different sound depending on the "phase" between the harmonics.

Following the same example, there is another way to see the phase between these two signals. Take for example the case in which $\phi_0 = 0$:

As you can see, whenever the fundamental wave crosses the $t$-axis, its 3rd harmonic also does. On the other hand, whenever the first one reaches a maximum (a minimum), its 3rd harmonic reaches a minimum (a maximum). This could lead to one thinking that they are in anti-phase. Again, this depends on what your application is and what you are looking for doing this analysis, as well as on what you define as "phase".

As a sidenote, there is a typical example that uses the phase between two waves of different frequencies (or something like that). It's called PLL (or phase-locked loop). You can check that out if you are interested in the topic.


To be precise the group delay of a linear phase FIR filter is $(N-1)/2$ samples, where $N$ is the filter length (i.e. the number of taps). The group delay is constant for all frequencies, because the filter has a linear phase, i.e. its impulse response is symmetrical (or asymmetric). A linear phase means that all frequency components of the input signal experience the same delay, i.e. there are no phase distortions. So for a frequency selective filter (e.g., a low pass filter), if the input signal is in the passband of the filter, the output signal is approximately equal to the input signal delayed by the group delay of the filter. Note that in general FIR filters do not have a linear phase response. In this case, the group delay is a function of frequency.

EDIT:

Consider an odd length FIR filter with even symmetry ("type I filter"). In this case the group delay is an integer number of samples $K=(N-1)/2$. The filter has a linear phase response (because the group delay is the negative derivative of the phase) and the filter's frequency response can be written as

$$H(\omega)=A(\omega)e^{-jK\omega}$$

where $A(\omega)$ is the real-valued amplitude function, and $\phi(\omega)=-K\omega$ is the (linear) phase. If the input signal's Fourier transform is $X(\omega)$, then the Fourier transform of the output signal is given by

$$Y(\omega)=H(\omega)X(\omega)=A(\omega)X(\omega)e^{-jK\omega}\tag{1}$$

Now assume that the frequency content of the input signal is in the passband of the filter, i.e. in the frequency region where $A(\omega)\approx 1$ holds. Then from (1) the output signal's spectrum can be written as

$$Y(\omega)\approx X(\omega)e^{-jK\omega}\tag{2}$$

and, from the shift property of the Fourier transform, the output signal is (approximately) given by

$$y(n)\approx x(n-K)$$

i.e. $y(n)$ is a delayed version of the input signal, where the delay is given by the group delay $K$.


Inverse $\mathcal Z$-transform of rational functions

Given the $\mathcal Z$-transform :
$$H(z) = \frac{\left(1+\beta z^{-1}\right)\left(1+\beta z\right)}{\left(1+\alpha z^{-1}\right)\left(1+\alpha z\right)}$$

The inverse $\mathcal Z$-transform can be found by the method of partial fraction expansion guided with the Region of Convergence ROC associated with $H(z)$

I will proceed with the most fundamental method and apply necesary algebraic manipulations to get the required standard form as follows:

\begin{align}
H(z) &= \frac {\left(1 + \beta z^{-1}\right)\left(1 + \beta z\right)}{\left(1 + \alpha z^{-1}\right)\left(1 + \alpha z\right)} &\scriptstyle{\text{convert to negative power of $z$}}\\
 &= \frac{\left(1 + \beta z^{-1}\right)\left(z^{-1} + \beta\right)}{\left(1 + \alpha z^{-1}\right)\left(z^{-1} + \alpha\right)} &\scriptstyle{\text{convert to standard form of } (1-d_k z^{-1})}\\
&= \frac{\beta}{\alpha} \frac{\left(1 + \beta z^{-1}\right)\left(1 + \frac{1}{\beta} z^{-1}\right)}{\left(1 + \alpha z^{-1}\right)\left(1 + \frac{1}{\alpha}z^{-1}\right)} &\scriptstyle{\text{expand the brackets}}\\
 &= \frac{\beta}{\alpha} \frac{\left(1 + \left(\beta + \frac{1}{\beta}\right) z^{-1} + z^{-2}\right) }{\left(1 + \left(\alpha + \frac{1}{\alpha}\right) z^{-1} + z^{-2}\right)} &\scriptstyle{\text{apply the long division}}\\
&= \frac{\beta}{\alpha} \left( 1 + \frac{ \left(\beta + \frac{1}{\beta} - \alpha - \frac{1}{\alpha}\right)z^{-1} }{\left(1 + \alpha z^{-1}\right)\left(1 + \frac{1}{\alpha}z^{-1}\right)} \right) &\scriptstyle{\text{go further...}}\\
&= \frac{\beta}{\alpha} \left( 1 + \frac{A}{1 + \alpha z^{-1}} + \frac{B}{1 + \frac{1}{\alpha}z^{-1}} \right) &\scriptstyle{\text{where $A$ and $B$ are found from the right side quotient of $H(z)$}}\\
\end{align}

Now assuming $A$ and $B$ are found as functions of $\alpha$ and $\beta$, we shall define all possible ROC of $H(z)$ to find the corresponding inverse transform $h[n]$.

Now If $0<|\alpha|<1<\frac{1}{|\alpha|}$ and ROC is such that $|z|<|\alpha|$ then: $h[n]$ will be unstable and anti-causal as
$$h[n] = \frac{\beta}{\alpha} \delta[n] + \frac{\beta}{\alpha} \left(- A\left(-\alpha\right)^n \text{u}[-n-1] - B\left(- \frac{1}{\alpha}\right)^n \text{u}[-n-1] \right)$$
Or else if $0<|\alpha|<1<\frac{1}{|\alpha|}$ and ROC is such that $|\alpha| <|z|<\frac{1}{|\alpha|}$ then: $h[n]$ will be stable and two-sided as
$$h[n] = \frac{\beta}{\alpha} \delta[n] + \frac{\beta}{\alpha} \left( A(-\alpha)^n \text{u}[n] - B\left(- \frac{1}{\alpha}\right)^n \text{u}[-n-1] \right)$$
Or else if $0<|\alpha|<1<\frac{1}{|\alpha|}$ and ROC is such that $|z|>\frac{1}{|\alpha|}$ then: $h[n]$ will be unstable and causal as
$$h[n] = \frac{\beta}{\alpha} \delta[n] + \frac{\beta}{\alpha} \left( A(-\alpha)^n \text{u}[n] + B\left(- \frac{1}{\alpha}\right)^n \text{u}[n] \right)$$
The cases for $0<\frac 1\alpha<1<\alpha$ are simply the same, merely interchanging $\alpha$ with $\frac 1\alpha$


Concept of finding forced response from difference equation?

considering your equation,

$$k\times u[n] - 0.25\times k\times u[n-2] = 2\times u[n] + u[n-1]$$

the left half of the equation will result in a series as shown below:

$$\begin{array}{cccccccc} [k & k & \frac{3}{4}k & \frac{3}{4}k & \frac{3}{4}k & \frac{3}{4}k  & \ldots &]\end{array}$$

the right hand of the equation is
$$\begin{array}{cccccccc} [2 & 3 & 3 & 3 & 3 & 3 & \ldots & ]\end{array}$$

for these two series to be equal,
$$\begin{array}{cccccccc} k = 2 \\ k=3\\ \frac{3}{4}k=3 \\ \frac{3}{4}k=3 \\ \vdots & \end{array}$$

It is not possible to equate these two sequences without removing first two points. In that case, $k = 4$. The sequences are equal or the equation holds for $k = 4$ and $n \ge 2$.

What we are getting here is total solution, so it will have the natural response and forced response. The output get stabilized and reach the stable value only after $n \ge 2$. This 'transients' before $n < 2$ are due to the natural response.


FWGN model implementation vs explaination (textbook pages attached)

I assume that the function Doppler_spectrum computes the desired power spectrum of the noise. Since the power spectrum of a filtered noise process gets multiplied with the squared magnitude of the filter's frequency response, the filter's frequency response needs to be square root of the desired power spectrum (assuming white noise at the filter's input).

In formulas, if a random process $X(t)$ with power spectrum $S_X(\omega)$ is filtered by an LTI system with frequency response $H(\omega)$, the power spectrum $S_Y(\omega)$ of the output signal $Y(t)$ is given by

$$S_Y(\omega)=S_X(\omega)|H(\omega)|^2$$


Pitch Detection: avoiding frequency doubling / halving

so Adamski, what you have here is what we call "the Octave Problem" and is well-known and oft discussed with people who work on the pitch detection problem.

first, some definitions.  given a quasi-periodic signal, $x(t)$, where

$$ x(t) \approx x(t+T) $$

in some local region where $t \approx t_0$, we normally define "pitch", measured in octaves as the base-2 logarithm of the local fundamental frequency, $f_0$, of a note relative to some reference frequency $f_R$.  "pitch" is always about an interval, which means it's expressed relative to something.

$$ p = p_R + \log_2 \left( \frac{f_0}{f_R} \right) $$

the local fundamental frequency is the reciprocal of the local period

$$ f_0 \triangleq \frac{1}{T}$$

so if $p = p_R$ it's because the fundamental frequency $f_0$ is equal to the reference frequency $f_R$.  sometimes we use A440 as the reference, most often i use MIDI note 0 as the reference and measure pitch in semitones which are $\frac{1}{12}$ octave.  in the MIDI convention, middle C is MIDI note 60 and A440 is MIDI note 69.  then the definition of pitch (relative to MIDI note 0) is

$$ \begin{align}
p & = 69 + 12 \log_2 \left( \frac{f_0}{440 \text{Hz}} \right) \\
& = 0 + 12 \log_2 \left( \frac{f_0}{8.1758 \text{Hz}} \right) \\
\end{align} $$

so when $f_0$ = 8.1758 Hz (5 octaves below middle C), the pitch is "0" and when $f_0$ is an octave above that, the pitch is "12".

alright, so the basic issue regarding the octave problem is this:

(1)  given that $x(t)$ is fully periodic with period $T$, that is

$$ x(t) = x(t+T)  \quad \quad \forall t $$

then the same is also true that $x(t)$ is periodic with period $2T$ or $3T$ or any integer multiple of $T$.  that is easy to show.  the same is nearly true for quasi-periodic $x(t)$, so there are many candidates of possible periods $T$ to choose from, which one do we choose?  usually the candidate period $T$ that is smaller than all the other candidate periods is the one we choose.  in other words, A440 (A just above middle C) sounds like A440 with period $\frac{1}{440}$ second, but mathematically it is also periodic with period $\frac{2}{440}=\frac{1}{220}$ second, but we don't hear it as A220 (the A just below middle C).

(2)  now here's the rub:  suppose we add just a teeny bit of a 220 Hz waveform to our A440.  but this 220 Hz waveform is, say, 60 dB lower in level than the A440.  what are we going to hear and, mathematically, what is the period of this waveform?

so, your assignment for the day is to think about that and, maybe tomorrow i'll come back and tell you maybe one or two "secrets" (they belong to me, so i think i can tell them, but lacking contracting with me to develop an alg for you, i won't be too explicit) about how to deal with this, but you first have to understand the nature of the Octave Problem which you get a hint of by answering for yourself the question in (2).

stay tuned.


Longer preambles make it easier for the demodulator to lock onto the signal well.  This is due to getting more preamble signal energy and because longer preambles make it easier to detect stuff like carrier offset.

Barker sequences are nice for all the reasons you state, but they are also very limited.  The longest Barker sequence is 13 bits and there is only one of them.  The longest Barker sequence length that has two sequences is four, so you'd be limiting yourself to a four bit long preamble, which is not a very useful preamble.

You may say to yourself, "I could just use the same 13-bit Barker sequence for both the I and Q."  If you do that, though, you have made a BPSK constellation, only you've rotated it by $\frac{\pi}{4}$.

If you want to use a QPSK preamble, then I would use pseudo-random sequences from linear feedback shift registers (LFSR's).  You can find the polynomials for maximal length LFSR's here.


Why use PSD to analyze a signal instead of its Fourier Transform?

The power spectral density (PSD) is a natural measure of the signal's power content with respect to frequency. A central part of non-parametric signal processing is to provide a "best" estimate  of the "true" PSD from knowing only one or some "realizations" with finite length. By taking into account the influence of stationary random processes, it should be noted that the Fourier transform is not defined for processes with infinite energy.
However, by looking at the second-order properties (autocorrelation), the true PSD can be well-defined, and the Fourier transform of  the realization  (on a finite horizon) can be used to provide an estimate of it.


How to apply low-pass filter to an I/Q signal to get baseband signal?

You did not specify the bandwidth of your modulated signal, but for discussion let us assume it is 50 MHz. (Meaning all the signal of interest is in 50 MHz centered at your carrier at 1 GHz).

You may be downconverting to a lower frequency IF or doing a direct conversion to baseband in which case you will need to do quadrature downconversion to not have image issues. I will assume this is what you are doing which means you will be multiplying your real signal at 1 GHz with a sine and cosine using a local oscillator at 1 GHz (all digitally since you are in a sampled system, so an NCO).
After you downconvert with this multiplication each I and Q will have your signal of interest at baseband plus an image at 2GHz (the sum of the two frequencies at the multiplier inputs and this is what you want to filter. ("I" representing the real component of your complex baseband signal, and Q representing the imaginary component of your baseband signal). You need to eliminate the image at 2 GHz, but given your signal is only 50 MHz wide (+/-25 MHz at baseband), you would ultimately filter down to that bandwidth to eliminate interference and reduce noise, but to "see" your signal you could simply use a moving average to implement a very simple low pass filter given this significant separation between your baseband signal and the image.  A moving average filter would for a Sinc function in frequency with its first null at 1/T where T is the length of your average in time. So for example, if you wanted to place the first null at 500 MHz, you would average over 20 samples:

$T=\frac{1}{500e6}= 2$ ns. With a sampling rate of 10GHz, this is equivalent to 20 samples. What is convenient about this choice is that it will also have a null directly at 2 GHz which is the image you are trying to filter out!  So simply performing a moving average over 20 samples at the output of your mixer will significantly supress your image.  See the frequency response below where 1 on the frequency axis is half your sampling rate or 5 GHz.  In the example I gave with 50 MHz BW, the passband droop is negligible and the image supression (at .4 of the plot) would likely exceed 30 dB (I am eyeballing that) over your BW.  Point is that is it a quick filter to confirm your downconversion without having to get into an explanation on more detailed filter designs (which you would ultimately likely want).

Here are more details on the down-conversion process in case it helps:
To demodulate this signal directly to baseband, you would multiply it by $e^{-j\pi f t}$ where f = 1 GHz.  This is implemented with a sine and cosine at 1 GHz, and for best noise performance if you were concerned with that, you would do a full complex down-conversion as : $$(I_1+jQ_1)(I_2-jQ_2) = (I_1I_2+Q_1Q_2)+j(I_2Q_1-I_1Q_2)$$
Where

$(I_1+jQ_1)$ represents your modulated signal at 1GHz and

$(I_2+jQ_2)$ is your local oscillator for performing the down-conversion.

Since your 1 GHz signal is likely a real signal, to implement the lowest noise approach I suggest, you would need to split in quadrature before sampling, or alternatively implement the Hilbert transform digitally.  This also would not require a filter (although recommended of course to eliminate noise and interference out of band) as this is a full complex multiplier and therefore provides image rejection. Otherwise if you are ok with the 3 dB noise penalty and can filter out the image, $(I_1+jQ_1)$ would be real so $Q_1 = 0$, and the downconversion simplifies to

$$(I_1)(I2−jQ2)=(I_1I_2)-j(I_1Q_2)$$

Once you do that multiplication, you should see your signal of interest centered on baseband as $$I_{out} + jQ_{out}$$

Where
$$I_{out} = I_1I_2$$
and
$$Q_{out} = -I_1Q_2$$

Describing the I output ($I_{out}$) (Q would be similar) to see the image that would need to be filtered:

$I_1= cos(\omega_ct+\phi)$ where $\omega_c$ is the carrier frequency of the received signal with some arbitrary phase offset $\phi$ .

$I_2 = cos(\omega_ct)$ is the local oscillator I signal, tuned to the carrier frequency.

Using the cosine product relationship: $2cos(\alpha)cos(\beta)= cos(\alpha-\beta)+cos(\alpha+\beta)$

$I_{out}= I_1I_2 = cos(\omega_ct+\phi)cos(\omega_ct)= \frac{1}{2}(cos(\phi)+cos(2\omega_ct+\phi))$

Showing the high frequency component to be filtered out (in your case 2 GHz) and the signal at baseband with the phase offset as received.

Following the same process for Q results in

$Q_{out} = -I_1Q_2 = \frac{1}{2}(sin(\phi)+ sin(2\omega+\phi))$


Fast way to calculate a few points in convolution?

I'll write names of time domain sequences in lower case and their frequency-domain counterparts in upper case, for example $X = \operatorname{DTFT}(x)$, where $\operatorname{DTFT}$ stands for discrete-time Fourier transform. I chose this transform, because it is unaffected by zero-padding of the time-domain sequences, allowing equations that mix frequency-domain counterparts of sequences of different lengths. With the upper-case frequency-domain symbols, time-domain convolution is denoted by multiplication, and a time-domain delay of $N$ samples is denoted by $z^{-N},$ where $z = e^{i\omega},$ with $\omega$ as continuous frequency. No frequency-domain thinking is required here. You can mentally translate the equations while reading to time-domain convolutions, delays, and additions. Addition in frequency domain is addition in time domain.

The frequency-domain counterparts $X$ of $Y$ of time-domain sequences $x$ and $x$ each of length $2N,$ with $N=500,$ can be written in partitioned form $X = X_0 + X_1z^{-N}$ and $Y = Y_0 + Y_1z^{-N},$ where $\{X_0, X_1, Y_0, Y_1\}$ are frequency-domain counterparts of time-domain partitions $\{x_0, x_1, y_0, y_1\}$ each of length $N.$ The partitions are non-overlapping sub-sequences of $x$ and $y$, stripped of delay which is then reintroduced in the equations. Convolution of $x$ and $y$ can be written in frequency domain as:

$$\begin{align}&XY\\
= &(X_0 + X_1z^{-N})(Y_0 + Y_1z^{-N})\\
= &(X_0 + X_1z^{-N})Y_0 + (X_0 + X_1z^{-N})Y_1z^{-N}\\
= &X_0Y_0 + X_1Y_0z^{-N} + X_0Y_1z^{-N} + X_1Y_1z^{-N}z^{-N}\\
= &X_0Y_0 + (X_1Y_0 + X_0Y_1)z^{-N} + X_1Y_1z^{-2N}\end{align}\tag{1}$$

Now all convolutions between sequences are convolutions between partitions of length $N$. Each result of such a convolution is of length $2N-1$. Denoting the frequency-domain counterpart of the first $N$ samples of the convolution result by $H()$ as in "head" and the frequency-domain counterpart of the delay-stripped sub-sequence formed by the remaining $N-1$ samples by $T()$ as in "tail":

$$\begin{align}= &H(X_0Y_0) + T(X_0Y_0)z^{-N}\\
+ &\big(H(X_1Y_0) + T(X_1Y_0)z^{-N}\big)z^{-N} + \big(H(X_0Y_1) + T(X_0Y_1)z^{-N}\big)z^{-N}\\
+ &\big(H(X_1Y_1) + T(X_1Y_1)z^{-N}\big)z^{-2N}\\
= &H(X_0Y_0)\\
+ &\big(T(X_0Y_0) + H(X_1Y_0) + H(X_0Y_1)\big)z^{-N}\\
+ &\big(T(X_1Y_0) + T(X_0Y_1) + H(X_1Y_1)\big)z^{-2N}\\
+ &T(X_1Y_1)z^{-3N}\end{align}\tag{2}$$

Now we have neatly split the full convolution result into length $N$ partitions. You were interested in the partition of length $N$ beginning at index $2N,$ or in frequency domain:

$$T(X_1Y_0) + T(X_0Y_1) + H(X_1Y_1)$$

The time-domain counterparts of these three terms can can be calculated by three FFT convolutions of length $2N$ compared to your single FFT convolution of length $4N.$ I don't know if that is of any use. You'd still be throwing away a total of $3N$ samples of convolution results, as originally.

But wait! The time-domain counterpart of $T(X_0Y_1) + H(X_1Y_1)$ can be calculated by a circular FFT convolution of length $2N$ between $x$ and $y_1.$ First, writing that as linear (not circular) convolution:

$$\begin{align}&XY_1\\
= &X_0Y_1 + X_1Y_1z^{-N}\\
= &H(X_0Y_1) + T(X_0Y_1)z^{-N} + H(X_1Y_1)z^{-N} + T(X_1Y_1)z^{-2N}\\
= &H(X_0Y_1) + \big(T(X_0Y_1) + H(X_1Y_1)\big)z^{-N} + T(X_1Y_1)z^{-2N}\end{align}$$

Wrapping that around to make it circular, $z^{-2N}$ wraps to $z^{0} = 1:$

$$H(X_0Y_1) + T(X_1Y_1) + \big(T(X_0Y_1) + H(X_1Y_1)\big)z^{-N}$$

The second half of that has what we wanted, $T(X_0Y_1) + H(X_1Y_1).$ We still need to calculate the time-domain counterpart of $T(X1, Y0)$ by another length $2N$ FFT convolution so we are down to two FFT convolutions of length $2N,$ with only a total of $2N$ samples of the results wasted.

We can also confirm that user Fat32's suggestion of using a length $3N$ FFT convolution works, by wrapping the latter part of Eq. 2 around whereby $z^{-3N}$ becomes $z^{0} = 1$:

$$\begin{align}&H(X_0Y_0) + T(X_1Y_1)\\
+ &\big(T(X_0Y_0) + H(X_1Y_0) + H(X_0Y_1)\big)z^{-N}\\
+ &\big(T(X_1Y_0) + T(X_0Y_1) + H(X_1Y_1)\big)z^{-2N}\end{align}$$

The middle line has what we wanted. So a single length $3N$ FFT convolution is also an option, with only a total of $2N$ result samples wasted.


Checking linearity/time-invariance etc. for recursive system

You made no mistake. Following the hint in Dilip Sarwate's comment, you will be able to show the following input-output relation of the system:

$$y[n]=x[n-1]-x[n-2]+x[n-3]-+\ldots=\sum_{k=-\infty}^{n-1}(-1)^{n-k-1}x[k]\tag{1}$$

Using this relation, showing linearity will be much easier.

Furthermore, note that any system described by a linear difference equation with constant coefficients is linear and time-invariant (if the initial conditions are zero!). And your system obviously belongs to this class.


If it is a sinusoidal signal, there will be peak (among the frequency bins) in the frequency spectrum corresponding to the tone's frequency.

Ratio of the magnitude of this peak to the sum of the magnitudes of all other bins (which are noise) correspond to Signal to Noise Ratio.

But when its a non sinusoidal signal (like the one in your plot) you have to consider the relevant band of the signal instead of a single peak. you can quantify it by specifying a frequency bin corresponding to the mojor frquency component in the signal plus some leakage (related to bandwidth) into the nearby bins.

Please take a look in the following matlab code.

N               = 8192; % FFT length
leak            = 50;
% considering a leakage of signal energy to 50 bins on either side of major freq component

fft_s           = fft(inptSignal,N); % analysing freq spectrum
abs_fft_s       = abs(fft_s);

[~,p]           = max(abs_fft_s(1:N/2));
% Finding the peak in the freq spectrum

sigpos          = [p-leak:p+leak N-p-leak:N-p+leak];
% finding the bins corresponding to the signal around the major peak
% including leakage

sig_pow         = sum(abs_fft_s(sigpos)); % signal power = sum of magnitudes of bins conrresponding to signal
abs_fft_s([sigpos]) = 0; % making all bins corresponding to signal zero:==> what ever that remains is noise
noise_pow       = sum(abs_fft_s); % sum of rest of componenents == noise power

SNR             = 10*log10(sig_pow/noise_pow);


You are partly right. However It's not the filter coefficinets $a$,$b$, but the filter impulse response $h[n]$ which you have to multiply with the exponential term $e^{j w_0 n}$ to shift its frequency response to a desired center of frequency $w_0$. Note that for FIR filters the coefficients are such that $a=1$ and $b[n]=h[n]$; i.e., the impulse response equals the coefficients $b[k]$, but that's not so with the IIR filters.

Assume that you have an FIR filter with coefficients $a=1$ and $b[k]$ so that its impulse response is: $$h[n] = b[n]$$ And let the frequency response of the filter be $H(e^{j \omega})$, then the frequency response of the new filter $h_+[n] = e^{j w_0 n} h[n]$ will be
$$ h_+[n] = e^{j w_0 n} h[n] \implies H_+(e^{j \omega}) = H(e^{j (\omega - \omega_0)}) $$

and similarly the frequency response of the new filter $h_-[n] = e^{-j w_0 n} h[n]$ will be
$$ h_-[n] = e^{-j \omega_0 n} h[n] \implies H_-(e^{j \omega}) = H_(e^{j( \omega + \omega_0)}) $$

Note if you want a real impulse response filter, say $h_0[n]$, then you should add those two sub shifted filters:
$$ h_0[n] = h_+[n] + h_-[n] \implies H_0(e^{j \omega}) = H_+(e^{j \omega}) + H_-(e^{j \omega}) = H(e^{j (\omega - \omega_0)}) + H(e^{j (\omega + \omega_0)}) $$

You can achieve the same result by the following multiplication then;
$$ h_0[n] = 2\cos(\omega_0 n) h[n] \implies H_0(e^{j \omega}) = H(e^{j (\omega - \omega_0)}) + H(e^{j (\omega + \omega_0)}) $$

MATLAB implmenetation would be as follows:

h = fir1(64, 0.1 );     % prototype lowpass filter
L = length(h);

w0 = 0.5*pi;            % bandpass filter new frequency
n = [0:L-1];            % discrete-time index vector

h0 = 2*cos(w0*n).*h;    % new bandpass filter

% Display prototype lowpass FIR filter
figure,subplot(2,1,1)
stem(n,h);title('Prototype Lowpass filter h[n]');
subplot(2,1,2)
plot(linspace(-1,1,1024), abs(fftshift(fft(h,1024))));
title('DTFT magnitude |H(e^{j \omega})| of lowpass filter h[n]');

% Display new bandpass FIR filter
figure,subplot(2,1,1)
stem(n,h0);title('Prototype Lowpass filter h_0[n]');
subplot(2,1,2)
plot(linspace(-1,1,1024), abs(fftshift(fft(h0,1024))));
title('DTFT magnitude |H_0(e^{j \omega})| of bandpass filter h_0[n]');


Features of a non-periodic signal for comparison

It is just an idea but how about a mean of the terms in a sliding window.

x = zeros(1,250);
c = randi(50,1,15);
x(c) = randi(500,1,15);
y = x;
d = randi(15,1,3);
y(c(d)) = 0;

figure;
stem(x,'LineWidth',1); hold on; stem(y,'LineWidth',1);
legend('Original Signal','Signal with Missing Samples');

N = 5;
L = numel(x)/N;
[m, m1] = deal(zeros(1,L));

for i = 1:L
    mx(1,i) = mean(x((i-1)*N+1:i*N));
    my(1,i) = mean(y((i-1)*N+1:i*N));
end

figure;
stem(mx,'LineWidth',1); hold on; stem(my,'LineWidth',1);
legend('Mean of x', 'Mean of y');

As you can see the plots are very similar now. You can also check their cross-correlation:

figure;
stem(xcorr(mx,my),'LineWidth',2);

These are of course just ideas to get you started. Now you can play with the window size and instead of doing $1:N$, $N+1:2N$ you can try $1:N$, $2:N+1$ , etc. and replace mean with other statistical parameters to see which one works the best for you.

Hope it helps!


What is causing these spikes in the output of the low pass filter?

I don't have MatLab with the DSP toolbox, but you're calling filter wrong.  Use:

y = filter(Hd, [1], x);

to set the feedback taps properly to $[1]$.

You should also inspect your filter design by using

freqz(Hd, [1], 512, 8000);

Although my hunch is that amplitude spikes in the output are directly caused by the filtering of the phase discontinuities you intentionally added when performing FM modulation.

Really, no modern equipment uses frequency switching to perform FSK.  FSK emitted from modern equipment is almost always continuous phase.  Are you really dealing with an ancient piece of Bell 103 modem hardware?

An example Octave script that generates the complex I/Q samples of continuous phase FSK is here:

https://dsp.stackexchange.com/a/53834/28112

if you need an example of how to generate CPFSK.


Is there an algorithm to compute the phase for a single frequecy?

The problem can be formulated as (nonlinear) least-squares problem:

$$F(\phi) = \frac{1}{2}\sum^{n}_{i=1}\left[ A \cdot \sin(\omega i+\phi) - f_{i}(\omega) \right]^{2}$$

where $F(\phi)$ is the objective function to minimize with respect to $\phi$.

The derivative is very simple:

$$F'(\phi)=\sum^{n}_{i=1} A\cdot \cos(\omega i+\phi)  \left[A \cdot \sin(\omega i+\phi) - f_{i}(\omega)\right]$$

The above objective function can be minimized iteratively using Gradient descent method (first order approximation), Newton method, Gauss-Newton method or Levenberg-Marquardt method (second order approximation - $F''(\phi)$ need to be provided in these).

Obviously, the above objective function has multiple minima because of periodicity, hence some penalty term can be added to discriminate other minima (for example, adding $\phi^{2}$ to the model equation). But I think the optimization will just converge to the nearest minima and you can update the result subtracting $2\pi k, k\in N  $.


[EDIT] After a second read, the proposed normalization looks non standard. Suppose that $m\le x \le M$ ($m$ and $M$ denote the min and max). The scaling factor will be, depending on the situation:

if $m\le M\le 0$: $-m$,
if $m\le 0\le M$, and $|m|\ge M$: $-m$
if $m\le 0\le M$, and $|m|\le M$: $M$
if $0\le m\le M$: $M$

It turns out to be (if I do not err) the largest of the absolute values of $m$ or $M$.
A more standard writing, more symmetric, could be:

x_max = max([abs(max(x)), abs(min(x))]);

or directly:

x_max = max(abs(x));

This computes the maximum x_max in absolute amplitude. This corresponds to an $\ell_\infty$ norm dispersion-type. Then you center $x$ around the mean (in some $\ell_2$ sense), then divide by the above maximum x_max.

If $\mu$ denotes the mean, you will end up:

for the first two cases, in $[\mu/m-1,\mu/m-M/m]$
for the last two cases, in $[m/M-\mu/M,1-\mu/M]$

And there is not a lot more I can say, without more information. If the extrema values are sufficiently symmetric around the mean, with a close-to-$0$ average, the outcomes would lay  around the $[-1,1]$ interval, because $M/m \approx -1 $ and $\mu/M$ or $\mu/m$ are close to $0$. But you could have different behaviors: if $m=-10$, $M=1$, $\mu=-9$, you will have x_norm in $[-0.1,1]$.
But if the mean is $0.5$, the interval would be $[-1.05,0.05]$.

Variable reduction  is something one could take classically as an input to a decorrelating or source separation technique;  a PCA for instance, or related tools (SVD, independent component analysis) for instance.
But subtracting the central location (here the average) before computing extremas and rescaling a dispersion after that is more standard.

And however, mixing different norms in normalization is not the practice I would recommend. So:

subtracting the mean and dividing by the (residual) standard deviation
subtracting the median and dividing by the MAD (median absolute deviation)
subtracting the midrange and dividing by the max or min of the residual (ie the half range)
mapping the values between fixed bounds, such as $[0,1]$

are much more common.


The polyphase components of $h[n]$ are

$$e_k[n]=h[Mn+k]\tag{1}$$

Note that if $h[n]$ starts at $n=0$, then so do all polyphase components $e_k[n]$. In order to get back the original sequence $h[n]$ from the polyphase components $e_k[n]$, we need to insert $M-1$ zeros between the samples of $e_k[n]$, and shift each component back to its original position, i.e. shift $e_k[n]$ by $k$ samples. In the frequency domain this corresponds to

$$H(z)=\sum_{k=0}^{M-1}E_k(z^M)z^{-k}\tag{2}$$

Filtering an input signal with $\mathcal{Z}$-transform $X(z)$ results in

$$Y(z)=H(z)X(z)=\sum_{k=0}^{M-1}X(z)E_k(z^M)z^{-k}\tag{3}$$

which corresponds to the structure in figure $4.37$ in Oppenheim and Schafer. The terms $E_k(z^M)z^{-k}$ can clearly be implemented either by delaying the input or the output of the polyphase components:

$$X(z)\left[E_k(z^M)z^{-k}\right]=\left[X(z)z^{-k}\right]E_k(z^M)$$


Explanation of binning in frequency analysis post-FFT

In its simplest form, the 'binning' process consists in summing the energies (squared magnitude) within groups of adjacent FFT values. This will give you the total energy in a set of disjoint frequency bands.

A more elaborate form consists in using overlapping triangle filter banks - you compute a weighted sum of the energy in a range of FFT bins to get a number which can be interpreted as the energy measured at the output of a band-pass filter of a given center frequency/width. For example, one can use this technique to build a visualization which uses a 1/3 octave scale on the frequency axis.


Wavelet image denoising: dual-tree versus double-density

The formalism of the continuous wavelet transform is relatively flexible. To make a practical tool out of it, you ought to discretize it, and here come the pain. It is quite easy to discretize it in a very redundant way (a discrete wavelet frame), but the price to pay when analyzing high-dimensional signals (2D images, 3D volumes) can quickly become too high: on the order or $J^d$, with $J$ the number of wavelet levels and $d$ the dimension.

So people have tried to obtain a critically sampled, discrete wavelet, without redundancy. This is the discrete wavelet transform (DWT). Nothing comes free. The price to pay is a narrow choice in the wavelets, and a processing sensitive to shifts. Ivan Selesnick, among others, worked on several intermediate designs: by allowing a bearable redundancy (typically $2\times$ or a little more in each dimension), one can obtain results close to, or even above the quality of fully redundant designs.

He worked especially:

on dual-tree wavelets and
on the double-density wavelet,

both  with a redundancy of $2$ in 1D, whatever the level. The first one compensates the poorer DWT choice by combining it with its Hilbert pair, allowing lower shift variance. The second allowed more degrees of freedom in the design of wavelet filters, hence a richer choice. Both have their merits, that may depend on the application. A little redundancy in wavelets often helps.

And both can be further combined, as in  The Double-Density Dual-Tree Discrete Wavelet Transform, by the same author in 2004. The gain comes from the double advantage  of having more design freedom in the wavelets (allowing for instance better cutoff frequency and symmetry) and the Hilbert pair for shift-invariance.

This is confirmed at Double-Density Complex Wavelet Transform:

  Both the double-density DWT and the dual-tree DWT have their own
  distinct characteristics and advantages, and as such, it was only
  natural to combine the two into one transform called the
  double-density complex (or double-density dual-tree) DWT.


For each single level of a isotropic or square DWT, the row/column order is not important.  Note however that there exists an alternative decomposition, termed anisotropic, rectangular, that does all rows first, then all columns (or the converse) which results in a different subband organization. This is described in the overview of 2D wavelets: A panorama on Multiscale Geometric Representations (arxiv) in Section 3: 2.3.2. Separable orthogonal wavelets.

[EDIT: Added references to related questions]

2D DWT Image Issue
What is the correct order of operations for a 2D Haar wavelet decomposition?
2D DWT computation order


My bitrate clock will have an instantaneous period of 32 or 40 ns (but a mean of 32.26ns)

That won't work great. I'd hardly call a 20% too long period "jitter", anymore. You're simply not synthesizing a useful clock; you're average correct rate doesn't help a receiver much to get a lock on your symbol timing.
The "easy" solution is use the facilities of your FPGA to actually produce a 31 MHz clock.
If there's no clock synthesizer that can do that for you: You'll have to resample your signal internally. There's no need that your sample rate needs to be a multiple of your symbol rate, it's just easy to deal with.
But really, you can calculate what your transmit signal would look like at any rate (e.g. sampled at 62 MS/s) and then resample it to a rate that your system actually works at (e.g. 62.5 MS/s). There's no need for an actual 15.5 MHz clock to generate a baseband signal at 15.5 Msym/s QPSK to transport your 31 Mbit/s.
So, what you need is a rational resampler with a rate of 125/31 or 250/31 or so.
Your FPGA vendor quite likely has an IP core for that. If not, it's no magic, but will take some development time. You're looking for a polyphase implementation of a rational resampler if you need to implement that yourself, so that you can first decimate your signal by 31 and then interpolate by 125 (or 250, or so).


How does modulation frequency appear in the modulation spectra?

This problem is precisely what synchrosqueezing wavelet transform was invented for, and indeed maps it with great precision. I'm still developing it, and first pre-release is expected today or tomorrow, but currently you can use it with all default settings to estimate the amplitude modulation information by plotting the appropriate row:

# printed output
1998  # carrier frequency
200   # modulator frequency

I'll be writing a detailed explanation of synchrosqueezing in context of CWT in a separate Q&A 'soon', but for now you can refer to sources in relevant functions' docstrings. The change to 200 (or less, relative to carrier_freq) was necessitated from limitations of working only with default values, but you can still use defaults with additional re-interpreting of max_row_fft.
As a tip, forget STFT for this task, CWT and CWT-based synchrosqueezing are superior. For more relevant intuition, see this answer (downvotes were due to a related 'fiasco').

How's it work?
I'll use STFT to explain; firstly you'll need more segments in total to capture the modulator; default segment length is 256, and overlap half of it, so to make more windows let's use noverlap=224 and don't touch nfft (no real need) - plot abs(Zxx), and the row at which it is maximal:

The modulator is seen more clearly now, and the plot on the right captures it (the AM "envelope"). If it's accurate, then the FFT of that row should produce a peak at the modulator's frequency - and it does:

To understand how this works one must understand first and foremost how the DFT works, and then the STFT; I'll leave good learning sources below, but also in a nutshell: STFT takes a windowed DFT over the signal, each window centered at a different point of the signal; if the window decays away from center, then signal points farther from center contribute less to spectrum's computation, thus localizing frequencies.
Each row of STFT uses the same window on a different frequency complex sinusoid (i.e. DFT basis), overall STFT forming a time-frequency representation of the signal. The row in STFT has an associated frequency, returned by scipy.signal.stft, which in this example is the carrier's frequency. As STFT values at that row capture the signal contents due to that row's associated frequency components, the row essentially is the original signal (in this case), thus taking its absolute value gives us the (carrier's amplitude) modulator signal, whose frequency we can find via DFT (FFT). One caveat, it's not just one row, as nearby frequencies also correlate to produce nonzero values; this is what synchrosqueezing addresses.

FT
DFT
STFT
DFT further intuition
DFT vs DTFT, zero-padding

"Row" clarification -- the "row" is the original signal by approximation (in this case - only one 'effective' row), except also having a 90-shifted imaginary component, so taking absolute value yields the enveloping (AM) signal. This is clear from a "continuous STFT" (shift window by 1 sample); it won't be exact (but note that the frequency is exact) as again, the rest of information is in other rows, and reconstruction involves extra steps (done by inverse STFT). -- Code

Code: (-- STFT version)
import numpy as np
import matplotlib.pyplot as plt
from ssqueezepy import synsq_cwt_fwd

#%%## AM signal generation ###################################################
duration = 1 # in seconds
fs = 44100 # Hz
carrier_freq = 2000 # Hz
mod_freq = 200
t = np.linspace(0,duration,endpoint=True,num=fs)
audio = (1 + 1 * np.cos(2*np.pi*mod_freq*t)) * np.cos(2*np.pi*carrier_freq*t)

#%%###########################################################################
Tx, fs, *_ = synsq_cwt_fwd(audio, t=t)
aTx = np.abs(Tx)

#%%## Plot synchrosqueezed CWT ###############################################
plt.imshow(np.flipud(aTx), aspect='auto', cmap='bone')
yt = ["%.2f" % f for f in fs[np.linspace(0, len(fs) - 1, 6).astype('int32')]][::-1]
plt.gca().set_yticklabels(['', *yt])
plt.show()

#%%## Find carrier freq ######################################################
max_row_idx = np.where(aTx == aTx.max())[0]
max_row = aTx[max_row_idx].squeeze()

# print peak's frequency
print(fs[max_row_idx])
#%%###########################################################################
# plot amplitude modulator
plt.plot(max_row[:800]); plt.show()

# find modulator frequency
max_row_fft = np.abs(np.fft.rfft(max_row))
plt.plot(max_row_fft); plt.show()

# peak at 200; exclude dc term
peak_bin = np.argmax(max_row_fft[1:]) + 1
print(peak_bin)

Note: amplitude values in above plot are off due to mistaken normalization used in MATLAB repository; this'll be fixed.


What's this discrete latch operation called?

I figured this out, it's called a Sample and Hold. The matlab page on it talks about how to implement one in matlab.

In summary, an interpolate function is used to create it. There are a couple of interpolate options, like "next", "previous" and "nearest". The "extrap" is needed to extrapolate the waveform for times 8 and 9. Continuing my example:

x     = [2,3,2,4]
t     = [1,3,4,7]
t_new = 1:9

x_new = interp1(t,x,t_new,"previous","extrap")

Produces

x_new =
     2   2   3   2   2   2   4   4   4


AutoCorrelation using FFT of power spectrum

In principle you are doing the right thing. Step 4 should produce a real result unless there is a coding error. Sometimes you add up some residual imaginary part due to numerical noise but, if any, that should be very small.

Here is an example

%% random signal of length n
n = 128;
x = rand(n,1);
% zero pad
x = [x; zeros(n,1)];
% fft
fx = fft(x);
% mag sqaured
fa = fx.*conj(fx);
% inverse FFT
a = ifft(fa);
% plot
plot(-n+(1:2*n),circshift(a,-n));
answer = {'no','yes'};
fprintf('Is a Real?: %s \n',answer{isreal(a)+1});


How to obtain impulse response from the differential equation of a system?

It looks like your transfer function is correct, but there's a small mistake in your partial fraction expansion:

$$H(s)=\frac{2}{(s+4)(s+2)}=\frac{1}{s+2}-\frac{1}{s+4}\tag{1}$$

The corresponding impulse response is

$$h(t)=(e^{-2t}-e^{-4t})u(t)\tag{2}$$

The response to $x(t)=te^{-2t}u(t)$ is indeed most easily computed by solving the convolution integral:

$$y(t)=u(t)\int_0^tx(\tau)h(t-\tau)d\tau\tag{3}$$

I leave the exercise of solving $(3)$ up to you, but if I'm not mistaken the result should be

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


[EDITED] From left to right, try to find the impulse response. On the bottom branch, delay the unit pulse $d$ and get $d[n-1]$. Add $Ad[n]$. Delay twice and get $d[n-3] +Ad[n-2]$. Add $Bd[n]$. Add $Cd[n-1]$. So the FIR coefficients of the impulse response $h$ will be $[1,A,C,B]$, but in the reverse order (as correctly pointed out by colleagues), because $y[n] = \sum d[k]h[n-k]$. So your system is causal, and $h[0]=B$, $h[1]=C$, $h[2]=A$, $h[3]=1$.

However, at $n=20$, you are in a stabilized region of the signal (all the known samples are equal to $1$, so even if you make a mistake as I did first, the answer will be the same, $1+A+C+B$, as stated before by other colleagues.


Triangulation from set of points from left camera

You don't need to estimate points in the second image, you need to match them.

One way is to take a patch of pixels around the point in one image and see which of the points in the second image has a patch most "alike".  One very simple way to do this is to take the sum of the squares of the difference between each pixel in the left-patch and the corresponding pixel in the right-patch.


2 point discrete differentiation is bound to produce highly noisy results.
try the 5-points stencil. you can also generate coefficients (i.e. more points) yourself using derivation of Lagrange polynomials.


Validity of an argument that two transfer functions are minimum-phase based on their ratio being minimum-phase

two transfer functions HL and HR can each be represented as a minimum-phase filter (MPF) plus a pure delay.

That is generally not true and it's easy enough to disprove it by counter example. Let's simply look at a first order all pass
$$H(z) = \frac{p-1/z}{1-p/z} = H_{MP}(z) \cdot H_{AP}(z)$$ which has a pole at $p$ and a zero at $1/p$. The "minimum phase" part, $H_{MP}(z)$ is simply unity, i.e. $H_{MP}(z)=1$ and the all pass part is just the original all pass. The all pass phase response is NOT that a of a pure  delay, we clearly have $\Phi_{MP}(\omega) \neq \omega \cdot t_0$
I think you misunderstand the point in the paper: I think they are just saying that a specific class of transfer functions can be approximated well with this model. In this case this class is "HRTF ratios for human heads" and they show that this works particularly well for spherical heads. Since binaural perception is mainly based on interaural differences, approximating the frequency dependent interaural time delay by a single constant value can be quite useful.


Count heart beat from audio capture using smartphone

This is going to be a long answer.

I'm going to start with an analysis of your files and what you can do to improve the audio.  I'll also try your files out on an algorithm I designed to detect heartbeats.  After that, I'll describe what went into designing my algorithm and how it works.  Then I'll make some general suggestions for how to get started with your own algorithm.

Analysis of existing files

I loaded each of your files into Audacity to listen to them and to see what the frequency content looks like.  In general, you don't merely have a problem with static (white) noise.  You have a lot of trouble with constant interference from what looks like the AC power frequency in your area (50Hz and harmonics.)

From worst to best:

File: static_with_beats_doesnt_evaluate.wav

Heavy interference from powerline noise.  There is a heartbeat signal in there that I can see.  It is too covered up with interference to be detected, though.  Nothing much will save this recording.  The 50Hz interference is too close to the actual signal to be separated.

File: in_accurate_with_noise.png
Here again your biggest problem is interference from the 50Hz power line.  At points, the signal is clean enough to pick out a heartbeat.  If I run it through a very sharp low pass filter at 45Hz, the algorithm I use can pick out individual beats but even then there are stretches where the heartbeat is simply inaudible.

File: accurate.wav

This file doesn't have the power line interference, probably because it was recorded outdoors.  It can be deciphered much better than the others, although there are still sections where there is simply no signal to be detected.  A steep low pass filter with a cutoff at 45 Hz gets almost perfect results from my algorithm.  There are short sections where the signal is too weak to pick out of the remaining noise.

File: accurate_72bpm.wav

This file is also missing the power line interference.  Occasional disturbances cause a beat to be lost here and there, but otherwise a very good recording.

Problems involved in detecting heartbeats

There are several basic problems involved in detecting heartbeats:

What is a heartbeat?  A heartbeat is composed of several sounds, most of which are hard to separate from background noises and interference.  You have to pick a "feature" than can be reliably picked out of the garbage.  Simply detecting the volume of the sound won't do (too many other things can be loud, too.)  What I've found is that the "lub dub" sound is a short burst of frequencies between about 30Hz to 40Hz.  That's what my algorithm looks for.
Threshold - When is have you really detected your feature that says there's been a heartbeat?  The volume changes continuously, and so does the volume of the noise and interference.  As fixed threshold won't do.  Your threshold has to adapt to the background and the changed volume of the heartbeats.  The volume of the heartbeats changes with each breath, and it changes as your microphone moves around on the chest.  I have a really squirrely solution for the threshold problem.
Noise - The algorithm must work in fairly wide range of noisy conditions.  The feature you look for should be fairly simple to pick out of nearly any noise.  The algorithm I use does well up to about 20dB signal to noise ratio.  It is insensitive to the absolute level of the heartbeat signal, but requires a fairly clean signal.
Hardware - A functioning algorithm does you no good if it depends on the sensitivity and frequency response of the microphone.  Whatever you do needs to work with any microphone on any phone.  The algorithm I use is very nearly completely insensitive to the hardware.  If the heartbeat is there, it will be found.

Point 3 is really the most important.  Loud doesn't matter, absolute level doesn't matter.  The whole thing lives and dies on how clean your signal is.  Regardless of what algorithm you use, a clean signal is a must.

No interference.  If you are picking up power line hum, fix it.  It's usually an electrical problem and can be solved with proper wiring and grounding.
Reduce motion noise.  In my experiments, I use an elastic band to hold the microphone on my chest.
Reduce background noise.  I use a small cardboard "biscuit" stuffed with cotton between the microphone and the elastic band.  It is made of several rings of cardboard with a plastic covering. The back is also cardboard.  The side facing my chest has only a layer of plastic between the cotton balls and the microphone.  That greatly reduces the background noise.

How I detect heartbeats

The method I eventually came up with is built around some observations of some (really) clean recordings of heartbeats.  What I found is that there's a sort of "burst" of frequencies for each "thump" of the heartbeat.  If you filter out everything but the content of the burst, then you get a fairly reliable heartbeat detector.  I use a bandpass filter for 24 to 44 Hz.  That seems to deliver the best results.  The range might depend on the person - I don't have enough samples to compare.

After the bandpass, I use a Hilbert transformation to get a continuous estimate of how loud the burst frequencies are.  The Hilbert function gives two outputs that are rotated 90 degrees from each other.  Square each, take the square root of the sum, and you have a continuous signal representing the amplitude of the "burst" frequencies.  Apply a low pass filter to that to get rid of wild variations.

Getting the correct threshold was the real problem, and I tried several ways.  I started with manually adjusting the threshold, but of course it changes as you breath and move so that didn't work well.  I tried an automatic calibration, which helped but had to be repeated anytime the microphone moved.  I tried all kinds of adaptive thresholds, and they all failed for various reasons.

What finally worked really well looks just absolutely stupid.

Take the log of the absolute value of the amplitude signal. Yes, take the logarithm of the amplitude (it doesn't matter what base logarithm you use.)  Pass that through a high pass filter, and apply a fixed threshold to the output.  That fixed threshold is not setting the absolute level of the signal.  It is setting the signal to noise ratio.  A lower value means "work with a lower signal to noise ratio" and a higher value means "require a higher signal to noise ratio."  The end effect is "trigger if there is a frequency between 24Hz and 44Hz present and it is some fixed ratio louder than the current noise level."

This setup is pretty much insensitive to the absolute levels of the input signal.  The frequency response of the microphone and the amplifiers pretty much doesn't matter.  It adapts automatically for movement when breathing and when the microphone changes position.

This diagram shows the central part of the process:

That's a snippet of a Pure Data patch that I put together to develop the algorithm.  The whole thing together with some example files is available on GitHub.  If you want to experiment with it, you will need Pure Data Extended 0.43 to use the patch - it depends on several external libraries that PD Extended includes but which other versions don't.

Summary

Regardless of how you choose to implement your algorithm, you have to start by looking at the available signal in your recordings.

Look at them in the time domain (oscilloscope view) and the frequency domain (FFT.)  Find out what your noise and your interference looks like.

Look at clean recordings and find characteristics that are different from the noise and interference.  Find something that can be picked out with a filter of some kind.

Design a filter that can separate the characteristic from the noise.

Devise a method for adapting the threshold to the current noise level and signal level.

My method has drawbacks, of course.  The biggest being that if there's no heartbeat it will "imagine" it is hearing one and start delivering false pulses.  The independence from the absolute level and hardware variations comes at the cost of having to guarantee a certain minimum signal to noise ratio.


How does one calculate a steering matrix for DOA given real data and only antenna positions

Figured it out reading https://arxiv.org/pdf/2008.12323.pdf
Basically, you construct the matrix for all directions. so 360x180 for all elevations and azumoths.


how do you compute the channel gain from path loss index in wireless communication?

You will have to properly define path loss index. As soon as you have done that, you'll have a formula that links path loss index and attenuation/gain.


Frequency of a noisy signal with fixed width pulses

Not an easy problem to solve...

You could always cross-correlate the signal with a replica square wave and then apply simple thresholding to find the pulse locations. Then just count the peaks in a given interval to determine frequency.

I think the square wave nature of the signal could make the signal quite difficult to deal with in the frequency domain because of all the harmonics that comprise it.


Is there any alternative basis for a fourier-like transform?

The most important reason for the "Fourier supremacy", as you called it, has been pointed out by Jazzmaniac in his comment. The complex exponential $e^{st}$ with $s=\sigma+j\omega$ is an eigenfunction of the convolution operator:

$$e^{st}*f(t)=\int_{-\infty}^{\infty}e^{s(t-\tau)}f(\tau)d\tau=
e^{st}\int_{-\infty}^{\infty}e^{-s\tau}f(\tau)d\tau=e^{st}F(s)$$

where the eigenvalue $F(s)$, if it exists, is the Laplace transform of $f(t)$, and for $\sigma=0$, i.e. for $s=j\omega$, $F(j\omega)$ is the Fourier transform of $f(t)$ (again, if it exists). Even though there are several important differences between the Fourier and the Laplace transform, they are deeply related, and for answering your question their differences are not important. By analogy, the same is true for the discrete-time Fourier transform and for the $\mathcal{Z}$-transform in the discrete domain. So when in the remainder of my answer I mention the Fourier transform, you might as well replace it by Laplace transform or $\mathcal{Z}$-transform, and everything will remain valid.

So why should we be concerned with eigenfunctions of the convolution operator? Because convolution describes the input-output relation of linear time-invariant (LTI) systems, and these systems are often very good approximations of practically important systems such as frequency selective filters, integrators, differentiators, etc.

Apart from being an invaluable tool for the analysis of LTI systems, the Fourier transform also has an important physical interpretation as the spectrum, i.e. the frequency content, of a signal. Furthermore, the FFT is a highly efficient algorithm for the computation of its discrete version for finite length signals.

So why don't we just say we're doing fine with the Fourier transform and its siblings? Because there are applications where we are looking for other properties of a transform. One example would be time localization in signal analysis, since the basis functions of the Fourier transform have no time localization at all because of their infinite extent. In such cases other transforms have become widely used, such as the short-time Fourier transform and the wavelet transform. For discrete-time signals, filter banks provide very general signal expansions. They can be designed to match the desired properties of the signal expansion.

Since you mentioned the saw-tooth function as a possible basis function I would like to point out to you the Slant transform (mentioned in this answer), which is used in image coding. Furthermore, there is the related Slantlet transform, which is a special type of discrete wavelet transform with piecewise linear basis functions. You can see the basis function in Fig. 5 of this paper.


Decimation and filtering in the frequency domain

I do a lot of decimation in the frequency domain. Little details are important.

I assume you already know the basic rules for fast convolution: the FFT length N is equal to the data blocksize L plus the length of the filter impulse response M minus 1. Each operation uses L samples of new data plus M-1 samples of data from the old block.

Ensure that the impulse response of your lowpass filter is shifted to the front of your time domain buffer AND properly windowed to M samples before you take the forward FFT to get the frequency domain representation of your filter. This keeps the result from wrapping around in the time domain when you take the inverse FFT. (Remember you're actually doing circular convolution when you want linear convolution.)

Kaiser is by far my favorite window because of its tuning knob. Use a large enough value to push the sidelobes down, as they will alias into your output. I typically construct a zero-phase brick wall in the frequency domain, take the inverse transform, apply the window and shift to the front of the buffer, then take the forward FFT.

Also make sure that the main lobe of your frequency response is still essentially zero past the Nyquist limit of your decimated output. I.e., don't make the Kaiser parameter too big.

If you're also doing frequency shifting by rotating the frequency bins, remember that you have to shift by a number of bins that corresponds to one data block size. For example, if L=M-1, then you can translate by any even number of bins.

Don't worry about FFT blocksizes that aren't powers of two. Just pick a convenient size that doesn't have any large prime factors and FFTW3 will perform well. I use L=M-1=3840 to decimate 192 kHz by 4:1 to 48 kHz with a block time of 20 ms. That's a FFT blocksize of L+M-1 = 7680 = 2^9 * 3 * 5.


As you said, in the fourier transform, the frequency of the wave depends on its position in the array. But, the amplitude of the wave is not always the real component, it is actually the magnitude of the complex number ($\sqrt{real^2+imaginary^2}$).

If your initial image is perfectly symmetric along the exact center, you can get real coefficients, but it's something which seldom appears in image analysis.

Edit: Here's more details... Couldn't find the picture I was refering to...

The discrete Fourier transform, the basis of every FFT algorithm, is defined as such:
$$X_k = \sum_{n=0}^{N-1}x_n\cdot e^{-i2\pi kn/N}$$

From Euler's formula, a complex exponential can be separated in a complex sum of cosinus and sinus:
$$e^{ix}=\cos(x)+i\sin(x)$$

If we now look back at the Fourier transform, it can be separated into two sums:
$$X_k = \sum_{n=0}^{N-1}x_n\cos(-2\pi kn/N) + i \sum_{n=0}^{N-1}x_n\sin(-2\pi kn/N)$$
As a sidenote here, if your image (or signal) is perfectly symemtric about it's center, the cosinus will resonate strongly with your image while the sinus will not output anything. It will result a completly real image. It's the other way around for a perfectly anti-symmetric signal where one half is the negative of the other.

Now knowing this, it's straightforward to compute the amplitude of the resulting FFT, using the formula I gave in the first paragraph.

To get the frequency of the image itself is less obvious. We would expect the result to look somewhat like this, with the zero-frequency component at the center. It is, however, not the case.

The "standard" FFT algorithm will put the zero-frequency component in a corner (I can never remember which one, top-left or bottom-left!), such that the array indice [0,0] corresponds to $k_x=k_y=0$. But wait, once you cross N/2, the pixels represents negative frequencies, counting backwards. A good reminder is Python's scipy.fftpack.fftfreq() function found fftfreq().

My personal tips to take a quick look at an image's spectrum is to shift it then display the logarithm of it. When it's not shifted, it's hard to see asymetries. When looking at the non-log version, it's usually just a big black square with a small fuzzy grey dot...


Compressed Sensing Mathematical Concept in Signal Processing

If the signal $x$ is already sparse in the original/time domain, then it does not need to be transformed, in other words the transform $\mathcal{J}$ can be taken to be the identity. So $x$ and $X=\mathcal{J}x$ are equal. What we want is to recover, from observed $y$, the  unknown $x$, via the recovery of $X$, since $x=\mathcal{J}^{-1}X$. In that case, $A$ is just the "random" sensing matrix.

The above paper does not seem too complicated, and I am not sure there are many "much simpler" explanation. You can try:

a YouTube primer: Introduction to Compressed Sensing
some links: Understanding Compressed Sensing
and perhaps from Emmanuel Candès: slides (Compressive Sensing - A 25 Minute Tour) and a tutorial (An Introduction To Compressive Sampling [A sensing/sampling paradigm that goes against the common knowledge in data acquisition])


Length of padding for multiple samples stored in a single MP3 file

The temporal masking effect is minimal. If your samples do not have extreme differences in loudness, one frame of silence should be enough. Depending on sampling rate one frame is >20ms, that should be plenty.


If you're talking about affecting the Perfect Reconstruction (PR) properties of a wavelet transform using wavelets with the PR property, then it won't affect it at all.  I.e., you can transform the zero-padded signal, then apply the inverse transform to the result, and you will get the original (zero-padded) signal again.  Remove the zeros from the reconstructed signal and you have the (actual) original signal.

As far as the DWT is concerned, the zero-padded signal IS the original signal.  So if you keep track of how many zeros you added before the decomposition, then you can just remove the same number after reconstruction.

Of course, modifying any of the coefficients before reconstruction will introduce artifacts into the result, but this isn't related to the fact that you zero-padded the input signal.  However, it can result in aliasing effects in the area which should otherwise be all zeros (after reconstruction).  But if, after reconstruction, you simply remove those values which weren't part of the "original signal" (the signal before you zero-padded it), then it won't make any difference (there could still be aliasing effects in the rest of the signal if you change the coefficients before reconstruction).

So if you're using zero-padding for the input signal, it doesn't affect the "nature" of the signal at all - it just means that, as far as the DWT is concerned, your original signal is the zero-padded signal and not your actual original signal.  Keep track of how many and where (beginning or end) your zeros are and you can just remove those values after reconstruction.


If $f(t)$ is a real non-negative function with Fourier transform

$$F(\omega)=\int_{-\infty}^{\infty}f(t)e^{-j\omega t}dt\tag{1}$$

the following must hold:

$$|F(\omega)|=\left|\int_{-\infty}^{\infty}f(t)e^{-j\omega t}dt\right|\le\int_{-\infty}^{\infty}|f(t)|dt=\int_{-\infty}^{\infty}f(t)dt=F(0)$$

since $f(t)=|f(t)|$.


I can't help with the fisheye lens issue, but assuming that effect isn't too large, the mapping between the above images is a perspective transformation.

We can transform a point on the plane $z=1$ into 3D space using a $3\times 3$ matrix.

$$\left( \begin{array}{c} X\\ Y\\  Z\\ \end{array} \right) =
\left( \begin{array}{c} s_x\ \ h_x\ \ t_x \\
                         h_y\ \ s_y\ \ t_y \\
                         p_x\ \ p_y\ \   1 \\ \end{array} \right)
\left( \begin{array}{c} x\\ y\\  1\\ \end{array} \right)$$

$s$ is scaling
$h$ is shearing
$t$ is translation
$p$ is perspective tilting (positive number means positive side tilts away)

We can project this object in 3D space back onto the plane $z=1$ simply by dividing the coordinates by $Z$. This maps every point in space to the intersection of the plane $z=1$ and the line connecting the point and the origin. This is basically the process of taking a photograph.

$$X_{proj} = X/Z$$
$$Y_{proj} = Y/Z$$

What we want to do is find the reverse transformation that takes us from the projected coordinates back to the original coordinates where the points only existed in the plane $z=1$.

It turns out that we can solve for all of the transformation parameters using 4 points (the four corners of a rectangle, for example) and solving the following system of equations (see here for the derivation and some example MATLAB code).

Making certain assumptions like the object not changing size or shape ($s_x=s_y=1$, $h_x=h_y=0$) can simplify these equations.