Confusion regarding a library of Digital Signal Controller

The function you are refering to only implements an order-1 filter. There is another function for order-2 filter.

If you want to implement an order-3 filter, you will have to split it in an order-1 filter and an order-2 filter. There are many ways to do it, but cascading the filters, as we would do in the analog world, is a good solution, albeit not optimal IIRC.

As to your original question, since the function refers to an order-1 IIR filter, an order-1 Elliptic or Butterworth, will have the same coefficients if you design them with the same cut-off frequency.

However, starting with order 2 or more, the coefficients will not be the same.


What's the difference between the Gabor-Morlet wavelet transform and the constant-Q transform?

Simply speaking both the const-Q-transform and the Gabor-Morlet wavelet-transform are just continuous wavelet transforms. Or, more precisely, approximations thereof, as there will always be discretization issues in real applications.

A property of wavelet transforms is that they have build in the constant Q-factor property, or in other words logarithmic scaling.
Gabor and Morlet are just two names of a particular wavelet function (complex exponentials with a gaussian window) which is used most commonly.
The CQ-transform just uses another basis function/wavelet and has a special name attached to it, probably to some historical reason.


Using ICA on EEG signals for feature extraction

How do I know I am not extracting noise (they are all very similar)?

You should first find or record a well understood set of signals with which you can test your feature extraction code to verify your setup is correct and that the extracted components are as expected.

Then you'll have higher confidence that you're not extracting noise due to errors in your code and you can focus on tuning your algorithms or reproducing other's results.


Proof of DTFT equal to DFT when signal is periodic?

The proof is rather simple if you consider the counter example and I will leave it to the OP to work out the details with these hints that are hopefully helpful.
Understanding the premise of the Fourier Series Expansion is a good starting point: Any single valued analytic function from $0$ to $T$ in time can be decomposed into an infinite series of sinusoids each with a frequency of $n/T$ where $n$ is any integer. Here is the key point:
The result of summing all those components will result in a waveform that matches the original function as expected from $0$ to $T$, but also since each of those components repeat themselves over that interval, then the original function will also repeat over every additional interval of $T$ in time! Consider that any other frequency other than these provided will not repeat as such, then they would disrupt the resulting repetition so therefore cannot exist.
When the beginning of a waveform is not continuous with the end of the waveform, a discontinuity results in the repeated pattern, leading to ringing at the edges in the Fourier reconstruction known as the Gibbs effect. This is consistent with the fundamental requirement that the waveform needs to be single valued and analytic function.


Checking if u[n] and v[n] are orthogonal

It's essentially asking you to show that two signals that are orthogonal in the baseband (their original form) remain orthogonal after amplitude modulation process. In other words, amplitude modulation preserves orthogonality.

Think about the question this way. You have two signals that are orthogonal, meaning that

$$ \sum_{m=-\infty}^{\infty}x[m]y[m] = 0 $$

This is a bit ugly, since we have to deal with infinities. However, your problem also says that they both have bandwidth $B$, which means that in frequency domain they will have finite width.

Let's define $\mathscr{B}$ as

$$ \mathscr{B} = \frac{B}{F_s}$$

which means that in your discrete Fourier transform signals $X[k]$ and $Y[k]$, the only non-zero samples exist for $-\frac{\mathscr{B}}{2} \le k \le \frac{\mathscr{B}}{2}$. Since DFT is an orthogonal transform, signals orthogonal in time domain will remain orthogonal in frequency domain, i.e.

$$ \sum_{l=-\infty}^{\infty}X[l]Y[l] =\sum_{l=-\mathscr{B}/2}^{\mathscr{B}/2}X[l]Y[l] = \sum_{m=-\infty}^{\infty}x[m]y[m] = 0 $$

We know that multiplying a signal by a cosine shifts is away from DC in frequency domain, and the fact that the modulating frequency $f_c$ is larger than $B$ means that this won't cause any problems with overlapping modulated bands (think about what would happen if this weren't the case).

I think this gives you enough hits to take it from here.


What's wrong with this implementation of my IIR filter?

If you take the Z transform on both sides of the difference equation:

$$Y(z)(1-0.695z^{-1}+0.666z^{-2})=X(z)(1.001z^{-1}-0.0010155z^{-2})$$

$$G(z) = \frac{Y(z)}{X(z)}=\frac{1.001z^{-1}-0.0010155z^{-2}}{1-0.695z^{-1}+0.666z^{-2}}=\frac{2.002z^{-1}-0.002031z^{-2}}{2-1.39z^{-1}+1.332z^{-2}}$$

It does not match your transfer function. I think you can work backwards to get the right difference equation now (note the x coefficients go on the numerator and the y's on the denominator). The system you implemented has a pole outside the unit circle so it's unstable.


How to filter this signal $x= \sin(2\pi 10 t)+\sin(2 \pi 40 t)+\sin(2 \pi 50 t)+\sin(2 \pi 60 t)$ using comb filters in time domain?

I guess you are trying to design a nortch filter at 50Hz. You can try it in this way.

Generate a signal in frequency domain filter_f=ones(4200,1).
Then you set the points corresponding to 50Hz and -50Hz to 0.
Transform the signal to time domain.filter=ifftshift(ifft(filter_f)).
window the signal. filter_w=filter.*hamming(4200).
Pickthe center part points of filter_w. And that is the final nortch filter.


determining FEC or parity bits from percent error

Is it possible to know with say a reed-Solomon technique how many paity bits I would need to add to each block to move in the neighborhoos of 100%?

Normally yes, if you know that you will never have more than $e$ bit errors in a block of length $N$, then you could design a RS code that would fix them all.  The problem in your case is that 30% bit errors is a LOT of errors.  Assuming the errors are randomly distributed, no RS code will be able to fix that many errors.  The number of parity symbols you would need would exceed $N$.

The analysis goes as follows- RS codes work on symbols, not bits.  All of the RS codes that I have seen have symbol widths of 6 to 9 bits, with the most common symbol width being 8 (CD's use an 8-bit wide RS code).  At 30% bit errors, well over half of the symbols will have errors.  You need two parity symbols for each error that you want to correct, so when over half of the symbols have errors it is impossible to correct them.

You may think "why not make the symbols narrower, reducing the chance that they will contain a bit error?"  The problem is that the block length, $N$, is $2^m-1$, where $m$ is the RS symbol width.  If you reduce the symbol width to 1, then your block length is 1 and you can't correct any data.  Thus, no RS code can correct your errors.

I would say that your best bet is a low-rate Turbo code, though turbo codes cannot guarantee that they will fix your errors.


Estimate Filter Coefficients from the Result of Linear Convolution with a Known Signal

The idea here is to build the problem in its Matrix Form.
We have the filter $ h $ which is represented by the matrix $ H $ to represent Linear Convolution operation:

$$ H = \begin{bmatrix}
{h}_{1} & 0 & 0 & \ldots &  & 0 \\
{h}_{2} & {h}_{1} & 0 & \ldots &  & 0 \\
\vdots &  & \ddots &  &  & \vdots \\
0 & \ldots & {h}_{m} & {h}_{m - 1} & \ldots & {h}_{1}
\end{bmatrix} $$

This means that the Linear Convolution is given by:

$$ \boldsymbol{y} = H \boldsymbol{x} $$

This is the equivalent of vY = conv(vX, vH, 'full') in MATLAB.
It is clear that for $ \boldsymbol{x} \in \mathbb{R}^{n} $ and $ \boldsymbol{h} \in \mathbb{R}^{m} $ (Namely the filter has $ m $ taps) then $ H \in \mathbb{R}^{ \left( n + m - 1 \right) \times n } $.

The above matrix $ H $ is Toeplitx Matrix.
You may be able to generate it using GenerateToeplitzConvMatrix() I added to my answer repository.

Bu we're after the coefficients of $ \boldsymbol{h} $. The nice thing is we can write the above as:

$$ \boldsymbol{y} = X \boldsymbol{h} $$

Where $ X \in \mathbb{R}^{\left( n + m - 1 \right) \times m} $ is still a Toeplitz Matrix since the convolution (When shape is full in MATLAB notation) is commutative.

So once we have the matrix $ X $ we can solve the above equation like any other linear system:

$$ \hat{ \boldsymbol{h} } = \arg \min_{\boldsymbol{h}} {\left\| X \boldsymbol{h} - \boldsymbol{y} \right\|}_{2}^{2} = {\left( {X}^{T} X \right)}^{-1} {X}^{T} \boldsymbol{y} $$

But since $ X $ has a special structure, being Toeplitz, there must be a way to exploit this, right?
Indeed there is, solving it using the Levinson Recursion. You'll be able to find a solver in the code repository of this question as well (Also see my Levinson Recursion Implementation). The idea is similar to Hilmar's solution.

But I had another idea. Toeplitz Matrix are a generalization of Circulant Matrices. If $ X $ was a Circulant Matrix we could use the highly efficient fft() function to solve the problem.

So the question is how to implement Circular Convolution using Linear Convolution.
Namely, we want to adjust the matrix $ X $ and the output $ \boldsymbol{y} $ so the will represent Circular Convolution between $ \boldsymbol{x} $ to $ \boldsymbol{h} $.

Let's say we know how build those (Let's call them $ \bar{X} $ and $ \bar{ \boldsymbol{y} } $) then we have:

$$ \bar{ \boldsymbol{y} } = \bar{X} \boldsymbol{h} $$

Yet since $ \bar{X} $ it diagonalizable by the DFT Matrix then (See Eigen Values of Circulant Matrix):

$$ \bar{X} = {F}^{H} \Omega F \Rightarrow \bar{ \boldsymbol{y} } = {F}^{H} \Omega F \boldsymbol{h} \Rightarrow \boldsymbol{h} = {F}^{H} {\Omega}^{-1} F \bar{\boldsymbol{y}} $$

Where $ F $ is the DFT Matrix and $ \Omega $ is a Diagonal Matrix of the Eigen Values of the matrix $ \bar{X} $. Those Eigen Values are basically the DFT Transform of $ \boldsymbol{x} $.

So basically we can solve this using the DFT (fft()). In MATLAB code it will be something in the lines of:

vH = ifft(fft(vYBar) ./ fft(vX), 'symmetric');

The full MATLAB code is available on my StackExchange Signal Processing Q45879 GitHub Repository.
It includes implementation of all 3 methods and one can see that each yields a perfect solution as long the system is well defined.

Resources

Circulant Preconditioned Toeplitz Least Squares Iterations.
Conjugate Gradient Methods for Toeplitz Systems.
Scientific Applications of Iterative Toeplitz Solvers.
Iterative Methods for Linear Systems with Matrix Structure.
Least Squares and Total Least Squares Methods in Image Restoration.
Fast Iterative Methods for Least Squares Estimations.
A Weakly Stable Algorithm for General Toeplitz Systems.
Algorithms for Toeplitz Matrices with Applications to Image Deblurring.
A Superfast Toeplitz Solver with Improved Numerical Stability.
Super Fast Solution of Real Positive Definite Toeplitz Systems.
Fast Solution of Toeplitz and Cauchy Like Least Squares Problems (Code in Giuseppe Rodriguez Page - TPLS).
A Super Fast and Stable Toeplitz Direct Solver via Randomization and Hierarchically Semi Separable (HSS) Methods (Code in Jianlin Xia Page - HSS Software Package).


Carrier frequency offset (CFO) estimation using cyclic prefix

To understand the equation in question consider how frequency offset $f_\mathrm{off}$ affects the received signal $r(t)$. For simplicity the authors assume the noise-free case and a frequency flat channel, i.e. the received signal before frequency offset is equal to the transmitted signal $s(t)$. Frequency offset is a frequency shift and can be expressed in time domain by
$$
r(t) = s(t) \mathrm{e}^{j2\pi f_\mathrm{off}t}.
$$
The received signal is now sampled at $t = n/F_\mathrm{s},\, n\in\mathbb{Z}$, where $F_\mathrm{s}$ is the sampling frequency. This yields
$$
r(n/F_\mathrm{s}) = s(n/F_\mathrm{s}) \mathrm{e}^{j2\pi n f_\mathrm{off}/F_\mathrm{s}}.
$$
In the following we write $r_n$ and $s_n$ for $r(n/F_\mathrm{s})$ and $s(n/F_\mathrm{s})$, respectively. The authors introduce the normalized frequency offset $\beta = f_\mathrm{off}/\Delta f$ but they don't define what $\Delta f$ is. From the context I understand that it is the subcarrier spacing $\Delta f = F_\mathrm{s}/N$, where $N$ is the number of subcarriers. Using this relationship we can rewrite the above equation as a function of $\beta$:
$$
r_n = s_n \exp\left(j2\pi \frac{\beta}{N}n \right).
$$
Next we consider a received cyclic prefix sample at $n$ (see above) and the corresponding received payload data sample $N$ samples apart at $n+N$:
$$
r_{n+N} = s_{n+N} \exp\left(j2\pi \frac{\beta}{N}(n+N) \right)
= s_{n} \exp\left(j2\pi\beta + j2\pi \frac{\beta}{N}n \right)
$$
In the above equation we have exploited that the cyclic extension is just a copy of the last $N_\mathrm{g}$ samples of an OFDM symbol and therefore $s_{n+N} = s_{n}$ (for all $s_{n}$ that are in the cyclic prefix). We can now see that the phase difference between the cyclic prefix and its corresponding data samples is $\arg(r_n^* r_{n+N})=2\pi\beta$ which is the answer to your first question. It is now easy to obtain an estimate $\hat\beta$ of $\beta$
$$
\hat\beta = \frac{1}{2\pi}\arg(r_n^* r_{n+N})
$$
It's just an estimation because in reality the received signal is corrupted by noise and channel distortion. The estimation can be improved by not only comparing one cyclic prefix and data sample but all $N_\mathrm{g}$ cyclic prefix samples with their respective copies. This is done by adding all $r_n^* r_{n+N}$ in the cyclic prefix and answers your second question about how the formula for $\beta$ is derived.

I have to say that the paper you're referring to is poorly written and that I probably  wouldn't have understood it neither if I hadn't had some prior knowledge about this topic. Thus I suggest you also read the original work about this topic van de Beek et al.:

J. J. van de Beek, M. Sandell, and P. O. Borjesson, “ML estimation of time and frequency offset in OFDM systems,” IEEE Trans. Signal Process., vol. 45, no. 7, pp. 1800–1805, Jul. 1997.


In the early days of analog television broadcast with CRT display receivers, it was found that instead of sending a whole picture frame information (576 PAL video lines) at once at a frame rate of 25 Hz, it was visually better to send half the lines as even field (0,2,4,...,574) and odd field (1,3,5,..,575) at 50 Hz. This is called interlaced video 576i for PAL and 480i for NTSC. If you send all the lines at once (and at half the rate of 25 Hz) then it's called progressive and denoted as 576p and 480p. These days they are used for 720 or 1080 etc line HD or FHD or 4k formats in the same line.

The reason for better visual characteristic is due to the fact that faster field rate yields smoother motion and yet reduced vertical resolution is compansated by several factors.


Adding a constant to a bandpass-signal in transmission

You're right, from a pure information point of view, the constant just wastes transmit power and doesn't add anything to the signal that you'd care about – the receiver will filter it out.

However, you've got to realize that all AM audio transmission schemes are positively ancient. So, receiver's abilities to demodulate a signal were very limited.

Simple (analog) AM receivers are nothing but a frequency selective circuit (e.g. a tuned antenna) and a rectifier, and a low-pass filter – done. If you're off a bit when tuning the circuit, it doesn't really matter.

Now, in your stereo AM (never seen that, wouldn't do that, seems strange and wasteful in spectrum) transmitter,  it's very important to tune to the right frequency for the first and second channel to be correctly received.

Now, that requires something a simple receiver could "lock" on. And that's exactly what the constant is: It adds a constant peak in the spectrum. If you have something like a PLL, you can lock on to that (make the PLL's bandwidth small, so that it's not affected by the AM).

Modern receivers could be much smarter – it's not very hard to guess the "center" of a double-sideband AM spectrum without a pilot tone.


System Frame Number in LTE Specifications

System Frame Number (SFN) is really just a counter. As its name, it tells us the index of system frame modulo 1024, the index is from 0 to 1023 (10 bits).

SFN does not help the physical synchronization (frequency carrier, DFT starting frame, pilots, ...) because it is the information obtained after the waveform synchronization. It is decoded in Master Information Block (MIB) in Physical Broadcast channel (PBCH) as 8 bits. The SFN of MIB changes every 4 frames (2 bits) thus its value is 8+2 = 10 bits.

Nonetheless, SFN does help physical layer protocol synchronization, i.e. synchronization in term of scheduling, between DL physical channels and UL physical channels (just from release 13 for cat-M or NB-IOT UE with long repetition), or synchronization between UL physical channels (for CQI reporting, SRS scheduling, ...). For more details you can take a look at specfication 36.213 all version, the scheduling of SRS, CQI, ... all depend on $n_f$ which is SFN, i.e. absolute SFN modulo 1024.


Generally I like to compose my sinusoids using this format:

x = cos(2*pi*f/fs*(0:num_samps-1))

Depending on the FFT routine you will use it will provide either the onesided or twosided result. Calculating the domain of the frequency axis is as follows. It should be noted that num_samps and NFFT may not be the same.

% Nfft - FFT size
% Fs - Sampling frequency in Hz
% oneside - flag indicating FFT results will be onesided or two sided

% Compute the frequency vector
if oneside == true
  f = fs * (0:Nfft/2).'./Nfft;
else
  f = fs * (-Nfft/2:Nfft/2-1).'./Nfft;
end

Have a look at this paper for a nice explanation on spectrum estimation.


How to find the actual size of an existing object from the pixels of an image?

I really am no expert in cameras (nor in OpenCV in which I'm pretty sure everything can be handled automatically et efficiently without doing calculations yourself). But without knowing the distance between the lens and your object, that will be difficult.

Actually, knowing the focal length of your camera and the distance from the nodal point to the image plane, you can deduce the distance from the nodal point to the object. Then just use Thales' theorem and you can do exactly what you want.

The problem of this method is that it relies on how well the focus is done on the camera. And because of the varying depth of field of your camera, the magnitude of the error will vary too. AFAIK, this method should work pretty well for short distances. For longer ones, that can be very difficult as the focal point is near infinity and the depth of field is really large.

As an alternative, you can put an object of known size next to the object which size you are measuring. Thales' theorem again and you get your data.


DFT of Transient vs. Analytical Solution

That's the case. Assuming you have a band-limited continuous-time signal $x_c(t)$, then the relation between its (analytically derived) CTFT, $X_c(j\Omega)$, and the (numerically computed) DTFT, $X_d(e^{j\omega})$, (or actually its DFT $X_d[k]$ ) of the samples, $x_d[n] = x_c(n T_s)$, of it taken at a period $T_s$ will be :

$$ X_d(e^{j\omega}) = \frac{1}{T_s} X_c( \frac{\omega}{T_s}) ~~~, |\omega| < \pi $$

In practice you are computing the N-point DFT $X[k]$ of $x_d[n]$ which is the uniform samples of $X_d(e^{j \omega})|_{\omega = \frac{2\pi}{N} k}$, for $k=0,1,...,N-1$. And they have the same magnitude relation too.

$$ X[k] =X_d(e^{j\frac{2\pi}{N} k}) = \frac{1}{T_s} X_c( \frac{\omega}{T_s}) ~~~, |\omega| < \pi $$

In other words, you will have a weight of
$$ K = \frac{1}{T_s}$$

between the computed DFT / DTFT and the theoretical CTFT of the signal.


Is there any trade-off in numerical precision or speed?

Yes. For delays that are integer multiples of the sampling period method 1 is far superior: it's computationally efficient, it's bit-exact, it's easy to implement and it's almost fool proof.
Method 2 is computationally expensive, you need to pick an FFT length (which is not trivial) it's subject to numerical noise and it incurs significant latency.  It implements a circular shift, not just a delay and so you need to manage overlap, frame boundaries, etc which is easy to get wrong.
Non-integer (or fractional) delays are completely different can of worms and quite complicated. Suggested reading: Splitting the unit delay


I don't quite understand why you feel "harmonics" are relevant to this discussion.

A 5 Mhz has sine wave has no harmonics. If the signal has harmonics, it is not a sine wave any more but a different signal (rectangular, triangle, etc.)
For ANY signal: you need to determine the highest frequency that's in the signal and then chose the Nyquist frequency to be higher. Otherwise you get either loss of information or aliasing
Most ADC do have an anti-aliasing low-pass filter build in: if you have harmonics higher than Nyquist, the filter will just chop them off.

Points 2 and 3 apply to ANY signal: noise, music, video, OFDM, and also harmonic signals like rectangle our saw-tooth waves. There really is nothing special about "harmonics".


power spectrum of passband of ofdm is it right shape or not?

Doesn't look right. a 64-FFT with only the center 52 carriers occupied should look pretty "rectangular".
I'd recommend not doing the upconversion. OFDM is designed to be an equivalent-baseband system, so you can definitely investigate whether your OFDM works directly in baseband. It's way more intuitive that way.


position of null carriers in N=256 subcarier OFDM

Positioning of null carriers in an OFDM signal is a protocol-dependent piece of information. Specifications will often call for some null carriers around the edge of the band (to allow for rolloff in the antialiasing filters in the front end). Another common place you'll find a null carrier is at the band center. This may be found in specifications that are targeted at low-cost (i.e. zero-IF) receiver technology. Zero-IF receivers' baseband DC offsets can cause large spurious responses in the tuned band center, hence the use of a null carrier there.

So, the short answer to your question is that it's unanswerable: it's dependent on the specifics of the system you're working with.


Use of I/Q representation for unmodulated Rx Signal

If you are "observing" the source, this implies there is some sort of information you are looking to get out of it, whether it be the total background noise, interference levels etc.
Do you find "value and use" in the RF signal's magnitude versus time? What about the RF signal's phase versus time? The IQ representation gives us both of these and it simply is another way to represent complex numbers:
$$Ae^{j\theta} = A\cos(\theta)+jA\sin(\theta) = I + jQ$$
One reason complex numbers (I/Q representation) for RF signals are useful is that for real signals, the positive and negative frequency spectrums are complex conjugate symmetric as complex numbers. This means we can represent any real RF signal using just the positive (or negative) spectrum as an analytic baseband signal, describing the signal independently from the carrier frequency.


Why Fourier transform and Stockwell-transform retain the absolute phase information of one signal?

It simply means that for a signal:
$$
x(t) = a(t) e^{j\phi(t)}
$$
that the phase information $\phi(t)$ is completely preserved in the transform
domain.

Time-frequency representations that use a quadratic form do not necessarily have this property:
$$
x(t+\tau)x^*(t-\tau) = a(t+\tau)a(t-\tau)e^{j(\phi(t+\tau) - \phi(t-\tau))}
$$
so if $\phi$ is of the form:
$$
\phi(t) = f(t) + C
$$
where $C$ is a constant, then $C$ is not preserved in $x(t+\tau)x^*(t-\tau)$.

This is important because sometimes you are interested in the absolute phase of a signal, not just its relative phase.


Non-coherent receiver for DBPSK: Implementation details

1) A typical approach would be to have a frame alignment word (and by "word" I mean bit sequence) at the beginning of every frame.  You could time synchronize by looking for that FAW through cross-correlations.  Barker codes are a good choice because they have good autocorrelation properties.

The problem with the approach is that the data, being random, can look like the FAW sometimes.  You can use the periodicity of the FAW to screen out those false positives.  For instance, if your frame is 1000 bits long, you would look for multiple hits spaced 1000 bits apart.  Once you have found them you have both time and frame synchronization.

2) There is not a direct relationship between carrier frequency and bandwidth.  There is a direct relationship between bit rate and bandwidth, but there are other factors involved.  The easiest way to determine your bandwidth is to take an FFT of your signal and measure it.

3) Your question is unclear.

4) You haven't provided enough information, and "optimal" can mean a lot of different things.  If you want a good answer to this you should start a new thread with more details.

5) It sounds like you are adding noise via AWGN, so you know what your SNR is.  With BPSK Eb/No is equal to the SNR.


Understanding the Bilateral Filter - Neighbors and Sigma

I would start with the many resources on this site:

Is the Bilateral Filter a Solution of Some Variational Method?
How to Validate Bilateral Filter Implementation?
Comparison Between Guided Filter (Edge Preserving Filter) and Gaussian Filter.
What Is the Bilateral Filter Category: LPF, HPF, BPF or BSF?
Understanding the Parameters of the Bilateral Filter.

Regarding your questions, Let's address them one by one.

Q: How to Determine the Number of Neighboring Pixels?

The classic neighborhood used by Bilateral Filter is along the axis. So it is basically determined by Radius parameter. In order to make sense, the radius parameter is aligned with the $ {\sigma}_{s} $ parameter (Spatial Standard Deviation). Some examples of the connection between them can be found in:

How to Get Rid of Ripples from a Gradient Image of a Smoothed Image?
Parameters of Gaussian Kernel in the Context of Image Convolution.
How to Calculate Gaussian Kernel for a Small Support Size?
Gaussian Blur In MATLAB and Connection to Image Resolution.

Q: Are All Adjacent Pixels the Neighboring Pixels?

Well, Theoretically you can chose any way to set the neighborhood set (Circle, The whole image, etc...). But as I wrote above, the most common option it a rectangle (Square in most cases).

Q: What Are $ {\sigma}_{s} $ and $ {\sigma}_{r} $ For?

The parameters determines the range of values in the range and spatial domain.
Let's assume our image values are in the range [0, 1] and the spatial indices are integers. Assume $ {\sigma}_{r} = 0.1 $, then pixels with values which are farther than $ \sim 3 {\sigma}_{r} $ will have effective zero weight since $ {e}^{ -\frac{ {0.3}^{2} }{ 2 \cdot {0.1}^{2} } } = 0.0111 $ before normalization.

Same logic for the Spatial parameter $ {\sigma}_{s} $. If we set it to $ {\sigma}_{s} = 4 $ pixels which are more than 12 indices apart will have negligible weight. This is why we limit the radius of the neighborhood and match it to the spatial parameter as adding pixels which have zero weight is waste of computing resources.


Given that you have limited number of samples, you can perform a Discrete Fourier Transform and represent your sound in frequency domain (frequencies and phases). After that if you want to apply a delay of D degrees you can simply use the known Fourier transform below:

$$ F\left \{ f(t-\theta) \right \} \rightarrow \hat{f}(w) \times e^{jw\theta} $$

Assuming that your sound is a number of sinusoids added together; to apply same phase delay to all frequency components is multiply a complex number $ e^{jw\theta} $ where $\theta$ is your desired phase delay.

On the other hand if you want to apply the same time delay to all your frequency components you should calculate the appropriate phase angle shift to be applied to each frequency component. Assuming $t_{delay}$ is the delay you want to apply, you can obtain $\theta_w$ to apply with
$$ \theta_w = w\times t_{delay} = 2\pi f\times t_{delay} $$

In which, $\theta_w $ should be calculated for each frequency bin of your results of your DFT. And the reasoning behind this can be simply explained as below. Assume a signal with sinusoid components.
$$ f(t) = sin(2\pi f_1t) + sin(2\pi f_2t) $$
If I want to apply the same phase angle shift to both signals (I don't know why anyone would want but anyways):
$$ g(t) = sin(2\pi f_1t + \theta) + sin(2\pi f_2t + \theta) $$
And to apply a time delay:
$$ f(t-t_{delay}) = sin(2\pi f_1(t-t_{delay})) + sin(2\pi f_2(t-t_{delay})) $$
$$ f(t-t_{delay}) = sin(2\pi f_1t-2\pi f_1t_{delay}) + sin(2\pi f_2t-2\pi f_2t_{delay}) $$
As you see the phase angle shift $\theta$ you should apply is different for each frequency component.

edit: If you have the proper MATLAB license you could try checking out the doc and source of delayseq function (it's in Phased Array System Toolbox).


Cepstrum argument is called quefreency, which is in fact a time domain. So for example if you are looking for the fundamental frequency then you are searching for a peak in a specific range. In your case that's $[0.002; 0.2] \;sec$, which corresponds to $[50; 500] \; Hz$., knowing that $f=\dfrac{1}{t}$.

So when you are searching for fundamental frequency, you can find peak in quefrency domain and then convert it to corresponding frequency. You can also use a different approach. We know that our cepstrum is symmetrical around center (similarly to Fourier Transform), and also we have our sampling frequency $f_s$. Let's say that you have $N$ cepstrum bins, therefore distance between each of them is: $\dfrac{f_s}{N}$. That's why you can find the corresponding fundamental frequency $f_0$ by dividing sampling frequency by the index of a peak.

For example our peak is at index number $10$, and sampling frequency is $1000 Hz$. Corresponding $f_0$ is $100 Hz$. If peak is at bin $50$ (more to the right, meaning lower frequency - remember quefrency is like time), then $f_0 = 20$ Hz. The subraction of 1 in denominator is cause by fact that MATLAB is indexing from 1 not from 0, like most of the programming languages.


How to find out amplitude of harmonic after diff type of FT?

Lets find the spectrum of the signals $x_c(t) = A \cos(\Omega_0 t)$ and $x_d[n] = x_c(nT) = A\cos(\Omega T n)$

I-First consider the CTFT (continuous-time Fourier transform) $X(\Omega)$ of $x_c(t)$:
$$x_c(t) = A \cos(\Omega_0 t) = A\frac{e^{j\Omega_0 t} +e^{-j\Omega_0 t} }{2} \longleftrightarrow X_c(\Omega) = A \pi \delta(\Omega-\Omega_0) + A \pi \delta(\Omega+\Omega_0) $$ The CTFT spectrum involves two impulses at $\pm \Omega_0$. Their magnitudes are infinite whereas their weights are $A \pi$

II-Second consider the DTFT (discrete-time Fourier transform) $X_d(e^{j \omega})$ of $x_d[n] = A\cos(\Omega_0 n T)$

Approach-1: Considering directly the signal $x_d[n] = A\cos(\Omega_0 n T)$ where $\Omega_0 T = \omega_0$ and by using the DTFT properties:

$$x_d[n] = A \cos(\Omega_0 T n) = A\frac{e^{j\omega_0 t} +e^{-j\omega_0 n} }{2} \longleftrightarrow X_d(e^{j \omega}) = A \pi \delta(\omega-\omega_0) + A \pi \delta(\omega+\omega_0 ) $$ The DTFT spectrum involves two impulses (in $-\pi < \omega < \pi$) at $\pm \omega_0 = \Omega_0 T$. Their magnitudes are infinite and their weights are $A \pi$

Approach-2: Considering the relation between DTFT $X_c(\Omega)$ and CTFT $X_d(e^{j\omega})$ due to sampling $$\boxed{ X_d(e^{j\omega}) = \frac{1}{T} X_c(\frac{\omega}{T})}$$

Given that the CTFT $X_c(\Omega) = A \pi [\delta(\Omega-\Omega_0) + \delta(\Omega+\Omega_0)]$ we conclude on the DTFT as:
$$X_d(e^{j\omega}) =  \frac{1}{T} A \pi [\delta(\frac{\omega}{T}-\Omega_0) + \delta(\frac{\omega}{T}+\Omega_0)]$$

Recalling the property of the impulse as:
$$ \delta(\frac{\omega}{T}-\Omega_0) = \delta( \frac{\omega- \Omega_0 T }{T} ) = T \delta( \omega- \Omega_0 T )=T \delta( \omega- \omega_0 )$$

Yields for DTFT $X_d(e^{j\omega})$:
$$X_d(e^{j\omega}) =  \frac{1}{T} A \pi [\delta(\frac{\omega}{T}-\Omega_0) + \delta(\frac{\omega}{T}+\Omega_0)] = A \pi [\delta(\omega-\omega_0) + \delta(\omega + \omega_0)]$$ The same result of approach-1 is observed: The DTFT spectrum involves two impulses (in $-\pi < \omega < \pi$) at $\pm \omega_0 = \Omega_0 T$. Their magnitudes are infinite and their weights are $A \pi$

III-Third consider the $N$-point DFT of the signal $x_d[n] = A \cos(\Omega T n) = A\cos(\omega_0 n)$

Approach-1: Considering directly the $N$-point DFT sum of the discrete-time signal $x[n] = A\cos(\omega_0 n)$ for $n=0$ to $n=N-1$.

(NOTE for the simple analysis we shall assume there's an integer $k_0$ such that $\Omega_0 T = \omega_0 = \frac{2 \pi}{N} k_0$ is satisfied; i.e., the frequency $\omega_0$ resides exactly on a DFT bin frequency for the index $k = k_0$, for $0 \leq k < N$.)

$$X[k] = \sum_{k=0}^{N-1} A \cos(\omega_0 n) e^{-j\frac{2 \pi}{N} k n } = \frac{A}{2} \sum_{k=0}^{N-1}[e^{j\omega_0 n}+e^{-j\omega_0 n}]  e^{-j\frac{2 \pi}{N} k n }  $$ Since we assume that $\omega_0 = \frac{2\pi}{N} k_0$ then we get:

$$X[k] = \frac{A}{2} \sum_{k=0}^{N-1}[e^{j\frac{2\pi}{N} k_0 n}+e^{-j\frac{2\pi}{N} k_0 n}]  e^{-j\frac{2 \pi}{N} k n }   = \frac{A}{2} \sum_{k=0}^{N-1}e^{j\frac{2\pi}{N} (k_0 -k) n} + \frac{A}{2} \sum_{k=0}^{N-1}e^{-j\frac{2\pi}{N} (k_0 +k) n}$$

We can evaluate the sums using the relation
$$\sum_{k=0}^{N-1}e^{j\frac{2\pi}{N} (k-k_0) n} =
\begin{cases}
N  &, \text{ for } k = k_0 \\
0 &, k \neq k_0 \\
\end{cases}$$
Which yields for DFT $X[k]$ as:

$$X[k] = \begin{cases}
0 &, \text{ for } 0 \leq k < k_0 \\
\frac{A}{2} N &, \text{ for } k = k_0 \\
0 &, \text{ for } k_0 < k < N-k_0 \\
\frac{A}{2} N &, \text{ for } k = N-k_0 \\
0 &, \text{ for } N-k_0 < k < N-1 \\
\end{cases}$$

Approach-2 considering that the finite length signal $x_v[n]$ is obtained from the infinite length signal $x_d[n] = A\cos(\omega_0 n)$ by multiplying it with a window $w[n]$ of length $N$. And considering its effect on the DTFT of the windowed signal $x_v[n] = x_d[n] w[n]$ and considering that DFT $X_v[k]$ is obtained as the uniform samples of the DTFT $X_v(e^{j\omega})$ of $x_v[n]$:

Since $x_v[n] = x_d[n] w[n]$ then;
$$X_v(e^{j\omega}) = \frac{1}{2\pi} X_d(e^{j\omega}) \star W(e^{j \omega})$$

We previously computed $X_d(e^{j \omega})$ as:
$$ X_d(e^{j \omega}) = A \pi \delta(\omega \pm \omega_0)$$

And it can be shown that the DTFT of the window is:
$$W(e^{j \omega})= e^{-j \omega ( \frac{N-1}{2} ) } \frac{ \sin(\omega N/2) }{\sin(\omega /2) }$$

And the DTFT $X_v(e^{j \omega})$ becomes:
$$ X_v(e^{j \omega}) = \frac{1}{2\pi} X_d(e^{j \omega}) \star W(e^{j \omega}) = \frac{A}{2} W(e^{j (\omega-\omega_0)}) + \frac{A}{2} W(e^{j (\omega+\omega_0)}) $$

Finally the DFT $X_v[k]$ of $x_v[n]$ will be the samples of this $V(e^{j \omega})$ at the frequencies $w_k = 2\pi k /N$. Under this special sampling we have the following:
$$X_v[k] = X_v(e^{j \omega})|_{w_k = \frac{2 \pi}{N} k} = \frac{A}{2} W(e^{j (\frac{2 \pi}{N} k-\omega_0)}) + \frac{A}{2} W(e^{j (\frac{2 \pi}{N} k + \omega_0)}) $$

and replacing again $\omega_0$ with $\frac{2 \pi}{N} k_0$ we have:
$$X_v[k] = \frac{A}{2} W(e^{j (\frac{2 \pi}{N} (k - k_0))}) + \frac{A}{2} W(e^{j (\frac{2 \pi}{N} (k + k_0)}) $$

Similar to approach-1, and considering $W(0) = N$ we have
$$
W[k-k_0] = \begin{cases}
0 &, \text{ for } 0 \leq k < k_0 \\
N &, \text{ for } k = k_0 \\
0 &, \text{ for } k_0 < k < N-k_0 \\
N &, \text{ for } k = N-k_0 \\
0 &, \text{ for } N-k_0 < k < N-1 \\
\end{cases}
$$

And therefore we conclude for DFT $X_v[k]$ as
$$
X_v[k] = \begin{cases}
0 &, \text{ for } 0 \leq k < k_0 \\
\frac{A}{2}N &, \text{ for } k = k_0 \\
0 &, \text{ for } k_0 < k < N-k_0 \\
\frac{A}{2}N &, \text{ for } k = N-k_0 \\
0 &, \text{ for } N-k_0 < k < N-1 \\
\end{cases}
$$


Power Spectral Density from Probability Density Function

If the $x[n]$ are i.i.d. random variables with mean $\mu$ and variance $\sigma^2$, then the power spectral density of this discrete-time random process (which is effectively white noise plus a possibly nonzero mean $\mu$) is what you have calculated.  The shape of the common pdf of the random variables is irrelevant in all this except insofar as the shape determines the mean and variance; for example. you would get the same PSD if the $x[n]$ were Gaussian random variables with mean $\mu$ and variance $\sigma^2$.


Clarification on eigenfunction property from Oppenheim's Discrete-time signal processing, 3rd ed

This statement can only be understood given its context. Right before the sentence you quoted, we have this statement:

It is no surprise that the condition for existence of the frequency
response is the same as the condition for dominance of the steady-state solution.
Indeed, a complex exponential that exists for all $n$ can be thought of as one that is applied at $n = -\infty$.

Consequently, starting out at $n=-\infty$, if the system is stable, all transients must have died out at "finite $n$" because at finite $n$ the moment that the exponential was applied is infinitely long ago.
We also depend on stability because otherwise the eigenvalue $H(e^{j\omega})$ wouldn't exist, at least not in the classical sense.


How to create a test to measure end to end audio delay

I recently ecountered this problem too (about a month ago). This is what I found..

Network delay measuring protocol (NDMP) - little useful free-ware program. I used the same program to measure other VoIP call delay.

You can download it from here: http://www.unforgettable.dk/42.zip
ReadMe file is included in .zip file.


[2D fourier transform]: Most people can't explain this

The scaling property of the Fourier transform tells you that if  $\mathcal{F}(x(t)) = X(f)$, then  $\mathcal{F}(x(at)) = \frac{1}{|a|}X(\frac{f}{a})$, with $a\neq 0$. So when you dilate in time, you contract in frequency, and vice-versa.

Now, for continuous 1D signals, the Fourier transform of a uniform box will be a cardinal sine. If you dilate the box, the cardinal sine will be contracted. This could be coarsely observed from the location of its zeros.

Here, you deal with a 2D discrete signal. So the result should be a discrete (or periodic) version of a cardinal sine. Its spread in one direction can be visually estimated by the horizontal and vertical black lines, since your 2D discrete Fourier transform is separable.

Now look at one horizontal line in your red boxes. The blackest dots that separate the horizontal gray/white bumps have the same location on both images. Now look at several lines. The black dots align vertically. So you see that the vertical lines are spaced in the same manner on the top and the bottom, which is exactly what you expect.

The result in the vertical direction is the controverse. Horizontal black lines are closer to each other in the bottom row.

What is sometimes misleading in Fourier plots is that some features in the frequency domain look orthogonal to those in the spatial domain. You can see more of this effect in An Intuitive Explanation of Fourier Theory.

If you want to gain additional conviction, crop a line from each spatial image, and an horizontal line from the spectrum, and superimpose the plots.


Generating a sound where end frequency is 1/2 of the start frequency

Think of frequency in this context more in terms of angular velocity. If that velocity is constant, nothing changes from the usual picture. However, if the frequency changes, it is more correct to compute the actual angle that is input into the trig functions as an anti-derivative of the frequency function.

Here the interpolated frequency is

$$f(t)=(1-\tfrac tT)\,f_0+\tfrac tT\,f_T,$$

so the angle at time $t$ is

$$\phi(t)=2\pi\,\int_0^t f(s)\,ds=2\pi\, [t(1-\tfrac{t}{2T})\,f_0+\tfrac{t^2}{2T}]=\pi t\,[(2-\tfrac tT)f_0+\tfrac tT\,f_T]$$

In contrast, your formula is

$$\phi_o(t)=2\pi\,t\,f(t)=2\pi t\,[(1-\tfrac tT)f_0+\tfrac tT\,f_T]$$

has a derivative or local angular velocity of

$$\phi_o'(t)=2\pi\,[(1-2\tfrac tT)f_0+2\tfrac tT f_T]$$

which has a zero at $t=\frac T2\,\frac{f_0}{f_0-f_T}$, which given your data is located at $t=0.43264s$.

Note that the correct formula also reaches the local angular frequency of zero, but at the double of that time interval, so outside $[0,T]$.

Note on implementation: The last form of $\phi(t)$ can be constructed as

 pi*(i*dt)*(f0+interpolate(i,f0,fT,N))

with N=samples and dt=T/N, provided that

interpolate(0,f0,fT,N)==f0 and interpolate(N,f0,fT,N)==fT


Correct way to implement Discrete and Stationary Wavelet Transform?

In the Haar wavelet case, I just used the following NumPy:

# assume x is passed in as a variable
haar_transform = zeros(N)
i = arange(N/2, dtype=int)
haar_tranform[i]      = x[2*i] + x[2*i + 1]
haar_transform[i+N/2] = x[2*i] - x[2*i + 1]

I believe the mathematical description includes some scaling by $1/\sqrt{2}$.

I've learned that while academic papers may describe something very complicated there's a simple way to compute it.


IFFT of Microphone Array Directivity Response

Your plot shows the magnitude of the array gain vs angle of arrival for the broadside beampattern. Before you can take meaning IFFT of this sequence, you need to follow the following steps.

The data in the plot is the magnitude of the array gain, which is stripped of all phase information. For IFFT to give desired result, we need to restore the phase information too. Fortunately, this can be easily accomplished by just changing the sign of the alternate lobes, assuming the underlying curve is of $\sin(Mx)/\sin(x)$ form.
The gain of a linear array looking at angle $\theta_0$ is proportional to $\dfrac{\sin\left( Ma(\cos(\theta)-\cos(\theta_0))\right) }{\sin\left( a(\cos(\theta)-\cos(\theta_0))\right) }$. So it makes more sense to take the IFFT of the array gain sequence when it is plotted against $\cos(\theta)$ as opposed to just $\theta$, or angle (as it is shown in your plot).

Also, assuming that the measurement has been made at the design frequency (i.e. sensor spacing $d = \lambda/2$), the result of the IFFT would provide the sequence of sensor weights. The number of FFT bins would only affect the zero-padding around the actual sensors.

I took the liberty of trying the above steps out on a set of data shown below. Note that it is very similar to your top plot.

The second subplot shown the result as we restore the sign by inverting the sign of the alternating lobes.

At this point, let theta and gain represent the x and y data in this plot (after sign correction).

Next, the data is resampled such that the x-axis is linearly spaced in $\cos(\theta)$ domain. Note that since $\theta$ ranges from 0 to 180, $\cos\theta$ would range from -1 to 1.

NFFT = 2^8;
xr = linspace(-1,1,NFFT);
yr = interp1(cosd(theta), gain, xr, 'linear', 'extrap');

The choice of NFFT is arbitrary, but it must be greater than the number of sensors.

The result of resampling is shown in the third subplot above. Notice how it looks more like a sinc function now.

Finally, we are ready for the IFFT.

wgt = fftshift(ifft(yr, NFFT)); % weight of sensors
% plot
figure;
plot((-NFFT/2:NFFT/2-1), abs(wgt), 'o-', 'linewidth', 2);
grid on;
xlim(20*[-1, 1]);
xlabel('Sensor location (\lambda/2)');
ylabel('Sensor weight');
title('IFFT result (after FFTSHIFT)');

Based on the above plot, it appears that there were 13 sensors in the experiment.


Force Linear Phase for a FIR Filter Synthesized Using Berchin's FDLS?

Remember that, as its name suggests, Berchin's frequency-domain least squares method (FDLS) is really just a least-squares fit of a desired frequency-domain response. Recall that for an arbitrary transfer function of the form:

$$
H(z) = \frac{\sum_{k=0}^{M} b_k z^{-k}}{\sum_{k=1}^{N} a_k z^{-k}}
$$

you can write the system's frequency response by simply letting $z=e^{j\omega}$:

$$
H(\omega) = \frac{\sum_{k=0}^{M} b_k e^{-j\omega k}}{\sum_{k=1}^{N} a_k e^{-j\omega k}}
$$

This sort of system can be implemented as an IIR filter with $b_k$ as the feedforward (numerator) coefficients and $a_k$ as the feedback (denominator) coefficients. The FDLS method simply requires you to define the desired (complex-valued) frequency response at a collection of frequencies $\{\omega_k\}$. Evaluating the above expression for each frequency $\omega_k$ yields a system of linear equations in the $b_k$ and $a_k$ variables. You can then solve this using least-squares methods to yield the best coefficients for the desired frequency response and filter order (in the LS sense, of course).

If you have the requirement for exact linear phase (i.e. better linear phase response than you might get by just applying the method itself and specifying a desired linear phase response), then you could reformulate the above a bit, if you're willing to accept an FIR-only filter structure.

Recall that linear phase in the frequency domain corresponds to symmetry in the filter's time-domain impulse response. For FIR filters, this impulse response is the same as the filter's coefficients $b_k$.
Let $a_k=0\ \forall\ k$ in the above expression, yielding:
$$
H(\omega) = \sum_{k=0}^{M} b_k e^{-j\omega k}
$$
The symmetry constraint corresponds to a loss in the number of degrees of freedom in the linear system; instead of each of the coefficients being independent variables, the linear-phase constraint requires that coefficients that are symmetric about the impulse response's midpoint to be equal.
Rewrite the above equation taking these constraints into account. For example, for a type I FIR filter (odd length, even symmetry), then you would have something like:
$$
H(\omega) = \sum_{k=0}^{\frac{M}{2}-1} b_k e^{-j\omega k} + \sum_{k=\frac{M}{2}+1}^{M} b_{M-k} e^{-j\omega k} + b_{\frac{M}{2}} e^{-j\omega \frac{M}{2}}
$$
where the filter has $M+1$ taps. You would then solve the system for the reduced set of coefficients, then reconstruct the overall filter impulse response using the symmetry constraint:
$$
h[k] =
\begin{cases}
b_k, && 0 \le k \le \frac{M}{2} \\
b_{M-k}, && \frac{M}{2}+1 \le k \le M
\end{cases}
$$
The result is an FIR filter that will have linear phase.


Convolution that results in an all-zero sequence

Expanding on my comment when this question was closed, suppose that
$$x[n]=e^{j2\pi n/3},~\forall n,\quad y[n]=\begin{cases}e^{j4\pi n/3}& n=0,1,2,\\0, &\text{otherwise.}\end{cases}$$

Then,
$\displaystyle \sum_{n=-\infty}^\infty x[k−n]y[n]= \sum_{n=0}^2 e^{j2\pi (k−n)/3}e^{j4\pi n/3}
= e^{j2\pi k/3}\sum_{n=0}^2 e^{j2\pi n/3}=0.$

More generally,   $\displaystyle x[n]=e^{j2\pi n/N},~\forall n,\quad y[n]=\begin{cases}e^{j2\pi (N-n)/N}& n=0,1,2,\ldots, N-1\\0, &\text{otherwise,}\end{cases}$
works as a solution for sequences taking on $N$ values.
Of course, nitpickers might argue that y takes on $N+1$ values, not $N$
as required by the OP, but in the case $N = 3$, there is an easy workaround.
Define
$$z[n] = \begin{cases}1, & n = 0,\\
-e^{j2\pi/3}, & n = 1,\\0, &\text{otherwise,}\end{cases}$$
which is clearly a three-valued sequence, and note that
$$\displaystyle \sum_{n=-\infty}^\infty x[k−n]z[n]=
\sum_{n=0}^1 e^{j2\pi (k-n)/3}z(n) = e^{j2\pi k/3}\cdot 1
+ e^{j2\pi (k-1)/3}\cdot (-e^{j2\pi/3})
= e^{j2\pi k/3} - e^{j2\pi k/3} = 0.$$ The real pezzonovante nitpickers
who wish to complain that the three values that $z$ takes on are
not the same three values that $x$ takes on are invited to work out
a more satisfactory solution for themselves.

It is tempting to extend $y$ periodically
to be nonzero for all $n$ so as to totally avoid the nitpick, but unfortunately,
the infinite sum $\sum_{n=-\infty}^\infty e^{j2\pi n/3}$ that is the
result of the convolution is not absolutely
convergent.  Of course, since both sequences are of period $3$, their
 convolution could also be defined as a periodic convolution,
which is the sum of just $3$ terms, and we are back in business.
More generally, note that the periodic
convolution of any two rows of the $N\times N$ DFT matrix is $0$,
and so we can extend this notion to sequences that take on $N$
distinct values for any $N \geq 3$.  For another viewpoint on this
approach, see the answer by @Hilmar.


Designing Rayleigh Fading Wireless Channel with particular properties


OFDM pilot allocation block type and comb type

In both cases, you need to interpolate the channel between the pilots you've got.
Both cases are typically suboptimal, since they'd only work perfectly for (a) actual block-fading (which is a convenient model, but doesn't look like reality) or (b) for a channel that is perfectly interpolatable from just a few points of observation in frequency (but that would imply you have designed an OFDM system with too many subcarrier, and that has other downsides).

You'd usually need both, because over the course of a longer data transmission, the channel might change, so there's need for regularly estimating more than just a few subcarriers. And because you cannot just guess the subcarrier state from a pilot some (up to) $\frac Q2$ subcarriers away, in general.

Modern systems tend to solve this by having something like (a) at the beginning of a frame, and then distributing pilots in a pseudorandom, or in a diagonal way throughout the frame, so that you get updates on the channel behaviour over time. Due to having updated your channel state info underway, you can space your pilot blocks further apart (i.e. increase $Q$), thus not reducing your data symbol rate by a lot.

LTE does that. There's a lot of resources you can find when searching for "LTE pilots", but linking to any of these here is a bit futile, since I don't know which background you're coming from.


I have worked with 4096-point Radix-2 based FFT algorithms, and although one might think that we need to adjust the accumulation after every stage by rightshifting the accumulation by 1 bit.

But that leads to over scaling of the magnitude and if the input $x[n]$ has low amplitude, then the accumulation at few stages might become 0 after this division by 2.

I found the scaling 0xAAA as most optimal performing one, which means that we are adjusting the accumulation at alternate stage rather than after very stage.

0xAAA : 101010101010 -> Accumulated register is right-shifted by the bit at bit-position in 0xAAA corresponding to the stage.

You can try with 0xFFF also, and actually see that the FFT result is getting higher error at those FFT Coefficients which are having low amplitude.

Hope this helps.


2-sided Regions of convergence for Z transforms

With one pole you have only two possible regions of convergence, either $|z|<r$  (where $r$ is the pole radius), which corresponds to a left-sided sequence, or the region $|z|>r$, corresponding to a right-sided sequence. I think the misunderstanding lies in the concept of 'the ring formed by the one pole', because there is no such ring.


what is the theory behind matlab butter function?

The MATLAB function that you refer to uses this knowledge to calculate the coefficients

Briefly, the function you refer to goes something like this:

The transfer functions for analog Butterworth filters of any order are already known analytically. It selects the normalized transfer function polynomial for the order you select.
The normalized polynomial is scaled appropriately according to the desired cutoff frequency, giving an $s$-domain transfer function that has the response that you want.
The $s$-domain prototype filter is then mapped to a discrete-time approximation using the bilinear transform. This is the most commonly-used method of mapping a continuous-time filter to a discrete-time filter.
The result of the bilinear transform is a transfer function in the $z$-domain, from which the filter coefficients are easily found by inspection.


Can compressed sensing be used instead of intepolation for missing values?

Yes, at least in the above case it is possible.
Though it might not be computationally as cheap as other methods such as least squares based curve fitting methods.
I do not think injecting NaN gonna help, instead  let's look at the ignored data as dimension reduction feature of CS.
Here, the measurement matrix is basically an identity matrix that some of its row (120:130) have been removed:
A = eye(256);
A(120:130) = [];

Just be aware, for the case of measurement matrix you posed in your question, it might violate coherence requirement.
If we apply the above matrix, we cut away samples of 120 to 130. Now assume the resultant as measurement vector $y$ and $A$ as measurement matrix. This is how CS recovery algorithms work, e.g. in case of Fourier sparse signals, we have an incomplete set of measurements $y$ ,the goal is to find sparsest set $a$ such that if we apply A matrix, yields something close enough to $y$.
N = 2^8;
t=linspace(0,16,N);
x=sin(2*pi*t)+sin(2*pi*3*sqrt(2)*t+pi/4);
x=x';
Phi = eye(256);
Phi(10:17,:)=[];
y = Phi * x;
Psi=dctmtx(N)';  % Transposed DCT matrix is like inverted DCT
recovered=SolveOMP(Phi*Psi,y); % any solver available!
CS_intepolate=idct(recovered);
plot(CS_intepolate,'*'); hold on; plot(x,'r')

(blue line is CS interpolated and red is original signal)
BTW, as far as I understood, you suggest to mix the remaining data. See it as system of equations, blending the current equations into new ones, does not add any new information or constraint to the system, so it is not helping.
Also I believe there is connection with CS and data completion. You might find it interesting also.


Different approaches for partial image derivation

Non-zero pixels after derivation provide you with two informations: a potential location (where the pixel is), and a potential strength (in magnitude) of an edge-pixel. Linear edge detectors are linear filters. A $3\times3$ kernel adds a smoothing effect in the orthogonal direction. It will favor edges with a certain spatial extend, over short-lengthed isolated pixels that do not qualify for real-edges. In other words, an isolated $[0,1,0]^T$, a duplet $[1,1,0]^T$ will be affected a smaller magnitude than three vertical ones $[1,1,1]^T$ after a $3\times3$ kernel.

In different words, the Prewitt  $3\times 3$ separable gradient/smoothing kernel can be written as a (Cartesian) product of two $3\times 1$ and $1\times 3$ kernels:

$$
\left(
\begin{matrix}
 1 & 0 & -1 \\
 1 & 0 & -1  \\
 1 & 0 & -1 \\
\end{matrix}
\right) =
\left(
\begin{matrix}
 1  \\
 1   \\
 1 \\
\end{matrix}
\right) \times \left(
\begin{matrix}
 1 & 0 & -1 \\
\end{matrix}
\right)
$$

This scheme can be generalized: any smoothing filter in one direction, combined with any derivative in the other (orthogonal) direction, can form an horizontal or vertical square/rectangular image edge detector.  You can find a source in Smoothed Differentiation Filters for Images, Meer and Weiss, 1992:

  For the two-dimensional case we restrict ourselves to square
  neighborhoods. In such neighborhoods it is always possible to define a
  separable two-dimensional orthonormal basis built by the Cartesian
  product of two identical one-dimensional bases.

and in Derivative-based Operations or the wiki page Image derivatives
.

If you have clear vertical edges, with same pixel amplitudes, little noise, the $1\times 3$ filter will do the job. Otherwise, good $3\times 3$ are almost always preferable. An image is really a 2D thing, and processing it through 1D tools often wastes information.


Analogous filter to Kalman that maximizes mode (as opposed to minimizing variance)

You can use a recursive Bayes filtering (RBF) approach to do maximum a posteriori (MAP) estimation. While MMSE estimators look at the conditional mean, MAP estimators look at the conditional mode. See these notes for some details. In order to do RBF, you should look into particle filters / sequential monte carlo methods. There are plenty of books on this topic as well as review papers.

Alternatively, if you assume that your Kalman filter has white Gaussian noise independent for the process and measurement, and an independent Gaussian initial state, you can prove that the MAP estimator coincides with the Kalman filter. [ A Kalman filter is a RBF for the nice Gaussian case. ]

A good paper to look at is S. Godsill, A. Doucet, M. West's “Maximum a posteriori sequence estimation using Monte Carlo particle ﬁlters,” Ann. Inst. Stat. Math., vol. 52, no. 1, pp. 82–96, 2001.

The actual article is behind a paywall, but certainly, many articles on related topics have been written (You can just google Maximum A Posteriori estimation with particle filtering or recursive bayes filtering for details).


Signal Processing using Fourier Transform

Hey it takes a lot of time in writing equations also I don't know to type in LATEX.so I am sending you image file for second one


Is Phase Locked Loop is essential for BPSK signal reception?

There are three different kinds of synchronization in a passband digital communications system:

Carrier synchronization: the receiver needs to know the exact frequency and phase of the carrier used by the transmitter.
Symbol synchronization: the receiver needs to know the optimum instants to sample the matched filter's output.
Frame synchronization: the receiver needs to know where each transmitted "character" (or word, or byte) starts.

So, if you're doing passband BPSK, you'll need to perform all three. Carrier synchronization is typically done with a PLL (or a Costas loop) and you can find the details in several answers on this website.

Many textbooks ignore these subjects. A textbook that does a great job of teaching these three aspects of a receiver is "Telecommunications Breakdown" by Johnson and Sethares, an early version of which is available for free from the authors' website.

Having said that, there is a way to avoid doing carrier synchronization, at the cost of energy efficiency, by upconverting the baseband BPSK signal using AM DSB-LC instead of the more common AM DSB-SC. So, you first design your baseband signal: $$s_{BB}(t) = \sum_k a_k p(t-kT_P),$$ where $a_k \in \lbrace -1, 1 \rbrace$, $T_p$ is the pulse interval, and $p(t)$ is the pulse shape. Then, the DSB-LC passband signal is $$s_{PB}(t) = (A_m + s_{BB}(t)) A_c \cos(2\pi f_c t),$$ where $A_m > -\min(s_{BB}(t))$ (so that $(A_m + s_{BB}(t))>0$) and $A_c$ is the carrier amplitude. Then, the baseband signal can easily be recovered by a simple envelope detector (identical to conventional analog modulation).


DFT (FFT) of Non Uniformly Sampled Signal

I am trying to perform an FFT of a non-uniformly sample signal. My input comes from Simulink and PLECS which uses variable-time solver.

Assuming that the first column in the data provided is your time variable ($t$) in units of seconds, the variable-time solver introduces a very small jitter.

(All code in Octave)

x = csvread("somefile.csv");
time_variable = x(3:end); %Skipping the first few samples as the solver settles.
plot(diff(time_variable)); %Ideally, this should return a flat plot at the sampling period, but that is not what we have here.

This results in:

Where you are sampling with at least a sampling period of something like $Ts \approx 1.68e-07$ extending to $Ts \approx 7.65e-07$.

The period of a 300Hz sinusoid that you are interested in is something like $0.00333 \ldots$ and the worst case scenario given above introduces an uncertainty of $\approx 0.022 \%$. That is per-cent.

Now, I am not sure which waveforms you are trying to obtain the spectrum of, but the second one is the one with the most discontinuities (the control signal?) and its spectrum looks like:

Fs = 1/7.6587e-07;
NFFT = 512;
FreqVec = Fs.*(0:(NFFT-1))./(NFFT-1);
semilogy(abs(fft(x(3:end,2),NFFT));
xlabel("Frequency (Hz)");ylabel("Amplitude");grid on;

Where as you can see, I have used the Sampling Frequency ($Fs$) of the worst case scenario, again, assuming that the first column contains the time variable in seconds.

Hope this helps.


If you define the SNR as the ratio of the signal power and the noise power in dB, you have

$$SNR_{dB}=10\log \left(\frac{P_s}{P_w}\right)\tag{1}$$

where $P_s$ is the power of the desired signal and $P_w$ is the noise poiwer. If the noise $w$ has a mean of zero, then $P_w=\sigma^2_w=1$. From (1) (with $P_w=1)$ you get the desired value of $P_s$ for a given value of $SNR_{dB}$:

$$P_s=10^{SNR_{dB}/10}$$

In order to normalize the signal $x$ such that it has the desired power $P_s$ you first need to know its power $P_x$. Dividing $x$ by $\sqrt{P_x}$ will give you a unity power signal, which can then be multiplied by $\sqrt{P_s}$ in order to obtain the desired SNR:

$$s=\sqrt{P_s}\frac{x}{\sqrt{P_x}}=10^{SNR_{dB}/20}\frac{x}{\sqrt{P_x}}$$


You can readily use the associativity of the convolution to explain this fact. Assume you send a signal $x(t)$ through a chain of convolutions $h_i(t)$:

$$ y(t) = (((x(t)*h_1(t))*h_2(t))*h_3(t)) $$

By associativity, this becomes

$$ y(t) = x(t) * ((h_1(t)*h_2(t))*h_3(t)) $$

I.e. you can see this operation as if sending your signal $x(t)$ through a single filter with impulse response $h(t)=h_1(t)*h_2(t)*h_3(t)$.


Do we need Channel estimation prior to (time domain DFE or linear) adaptative equalization in single tone transmission?

Linear Equalizer: Depdends on the type, but no, generally you just estimate an optimal inverse ("optimal" according to the specific metric of linear equalizer type) of the channel impulse response and apply that.
Decision-Feedback Equalizer: No. The way the channel equalization is computed is literally in the name.


Frequency shifting of a quadrature mixed signal

If you have a quadrature mixed signal, you have two outputs of the signal in 90 degree phase to each other. To complete a full complex frequency conversion (which would shift the frequency in one direction, as opposed to a real frequency conversion which would create two sidebands) you need to multiply one of the outputs by the cosine term of the required frequency shift and the other output by the sine of the required frequency shift. By swapping sine and cosine ports, you change the direction of the shift.  The outputs after each multiplier are then summed in phase. See diagram and derivation using trig identities for the multiplication of sines and cosines below.

If you are comfortable with the idea of exponential frequencies (full complex frequency including positive and negative frequency etc) then I prefer the view point below to describe single sideband frequency translation. Comparing the two diagrams shows the value of converting "I" and "Q" signals as typically constructed in hardware to implement complex signals into one complex signal path for the purpose of analysis. Using cosines and sines instead of exp can be quite cumbersome especially as the complexity of the system increases.

For those less familiar with the notation, $Ke^{j\theta}$ is simply a complex value with magnitude K and angle $\theta$. Therefore $e^{j\omega t}$ is a signal whose magnitude is 1 and angle is increasing linearly with time at rate $\omega$; a spinning phasor counterclockwise on a complex IQ (polar) plot, representing a "positive" frequency. Similarly $e^{-j\omega t}$ spins clockwise and represents a "negative" frequency. With that, one should be well equipped to use the more compact and analytically simpler exponential notation.

Notice too how the two block diagrams that I included are related, with the second one taking the real of the product to make the outputs equivalent. To do frequency translation, we use complex conjugate multiplication (the conjugate assures we shift in the direction desired). Consider a full complex conjugate multiplication, as is done prior to the real operation:

$ e^{j\omega_ct} = I_1+jQ_1$     representing the carrier after the quadrature split

$ e^{-j\omega_st} = I_2-jQ_2$     representing the conjugate of the I and Q terms of your signal

$ e^{j\omega_ct}e^{-j\omega_st} = (I_1+jQ_1)(I_2-jQ_2) = (I_1I_2+Q_1Q_2)+j(I_2Q_1-I_1Q_2)$

The real portion is also called the dot product and the imaginary portion is called the cross product. Notice the dot product does match the implementation block diagram shown at the top of this answer.

Further details for the very interested to understand complex frequency translation:

The following additional figures and content were added to help answer the question on the subsequent down-conversion and to build a better understanding of the complex (positive and negative) frequency domain in describing frequency translation implementations.

First consider a generic quadrature mixed signal, it is at "baseband", meaning not modulated to a carrier for purpose of transmission, and centered about DC (0 Hz) with the positive and negative frequencies independent of each other. The way we implement such a signal in practice is with two real signals (we could choose the real signals to have one represent the real and the other the imaginary such as I and Q in I+jQ, or one represent the magnitude K and the other represent the phase $\theta$ such as $Ke^{j\theta}$. In any event, by observing the spectrum and noting that the positive and negative frequencies do not match, we know from that detail it is a complex signal and two real signal paths are required to represent this in implementation. For a real signal in contrast, the positive and negative frequencies would be conjugate symmetric, meaning identical in magnitude and opposite in phase.

Frequency spectrum of baseband quadrature modulated signal (this could be a single tone on one side, or multiple independent frequencies, point is the positive and negative frequencies are independent):

Now consider a carrier signal $cos(\omega_ct)$, this has an impulse in the frequency domain at the positive frequency $\omega_c$, and another at the negative frequency $-\omega_c$ (also shown as two complex frequencies using Euler's identity $2cos(\omega_ct) = e^{j.\omega_ct}+e^{-j\omega_ct}$):

If we split the carrier in quadrature (with a Hilbert transformer, or a 90° splitter, or just generate a sine and cosine component, but in the end have two tones with a 90° relationship to each other), then in implementation we have two real signals used to represent a single complex frequency $e^{j\omega_c}$  (you may start to see how it is so much easier doing all analysis with the e's and then do the implementation with the sines and cosines):

UPCONVERSION

To up-convert the baseband signal to the carrier frequency, we multiply the two signals in the time domain (which is convolution in the frequency domain; so for an impulse function such as our carrier it is a simple shift). First consider what would happen if we multiplied the baseband signal with the carrier $cos(\omega_ct)$ directly (eliminate the 90° splitter in my first figure, replacing $sin(\omega_ct)$ with $cos(\omega_ct)$). Notice the output would still be complex as the positive and negative frequencies are not the same- so we could not transmit this with a single antenna (2 real signals are required to represent a complex signal). If we proceeded to take the real portion,
distortion of our signal would result as depicted in the figure, so we can't transmit that either.

Figure showing spectrum result with an INCORRECT upconversion approach:

This explains why we need to also take the Hilbert transform of the carrier (converting the cosine into sine and cosine components, and representing a  complex version of the carrier that is a single spinning phasor $e^{j\omega_ct}$, or in the frequency domain a single impulse in the positive domain with no negative component) prior to multiplication. We take the real portion of the product and can then send an undistorted real waveform representing our baseband signal modulated up to the carrier frequency.

Figure showing spectrum result with a CORRECT upconversion approach:

DOWN-CONVERSION

With the details of the up-conversion process in mind, it will now be very easy to answer the question if a quadrature conversion of the received signal is required prior to multiplying by sine and cosine to recover the baseband signal. The answer is no, but there is an advantage to doing so in situations when an image frequency cannot be easily filtered.

Figure showing spectrum doing a downconversion with no Hilbert transform of the received signal. It is interesting to note if you swapped sine and cosine in the (real) multipliers, this would represent a single impulse in the positive frequency domain, shifting the negative frequency component of the received signal to baseband, causing the spectrum at baseband to be inverted:

Figure showing spectrum doing a downconversion with a Hilbert transform of the received signal prior to sine/cosine multiplication. Here you can see what would happen if you swapped an input to one of the (real) multipliers: The signal would be shifted to a higher frequency twice the carrier instead of to baseband.


How to filter a signal knowing the frequency response of a system in MATLAB?

Your signal is given by

$$x(t)=\sum_{1\le|k|\le 6}a_ke^{jk\omega_0t}\tag{1}$$

with $\omega_0=2\pi/T$.

An LTI system's response to $a_ke^{jk\omega_0t}$ simply is $H(jk\omega_0)a_ke^{jk\omega_0t}$, where $H(j\omega)$ is the system's frequency response. Consequently, the response to $x(t)$ is

$$y(t)=\sum_{1\le|k|\le 6}H(jk\omega_0)a_ke^{jk\omega_0t}\tag{2}$$

If you want to use Matlab, you could better use it to evaluate Eq. $(2)$, rather than to simulate the continuous-time system.


The bandwidth of the message content of a CW signal in Hz is roughly at least 4 times the WPM at which the dots and dashes of each Morse Code character are being keyed.  The actual signal may take up a several multiples of that bandwidth, depending on the envelope shaping or the rise/fall time of the CW transmitter used to prevent key clicks. (Plus a few harmonics of a Morse Code dot's AM modulation envelope may be required for a typical radio operator to reliably copy dots by ear, requiring a bandwidth multiple corresponding to that harmonic number).  So any audio filter bandwidth needs to be at least that wide to pass a Morse Code signal, either by measuring the WPM of the signal, or using some maximum WPM expected.


There's no known name for this kind of filter in the signal processing world, as far as I know – and what good would such a name be, anyway?

So, I'd describe it in a signal processing context as "discrete state moving window sum with threshold", since it's a non-linear filter (and not much of a signal filter, if you ask me!). Simply because that's how I'd implement it: Sum up the last $W$ states, and output a 1 iff the result is $W$.

If you leave the realm of signal processing and think more of the basic electronic circuits that implement that, or by the purpose: You could call it a watchdog, or a digital anti-ringing, or a self-resetting cycle counter with a threshold, for example.

Your approach with the latch is from the latter category of names; "latch" is not anything that I'd call a filter in a signal processing context. Especially, the functionality you describe is not a latch (alone).


Number of zeros of a sum of Shah functions by applying Rice's formula?

This is not a complete answer in that it does not apply Rice's formula although it does derive the formula for the counting function $\mathcal{Z}(T_{max})$.

For $t,t_0,T,T_{max}\in \mathbb{N}$:

$$\frac{1}{T}\sum_{n=1}^{T}\exp\left(\frac{i 2\pi n t}{T}\right)=\cases{1&$T|t$\cr 0 &otherwise\cr}$$
so:

$$S(t,T_{\max})=\sum_{T=1}^{T_{\max}}\frac{1}{T}\sum_{n=1}^{T}\exp\left(\frac{i 2\pi n t}{T}\right)$$
counts the number of divisors of $t$ that are $\le T_{max}$ and $t_0$ is defined to be a zero if it satisfies:
$$\sum_{T=1}^{T_{\max}}\frac{1}{T}\sum_{n=1}^{T}\exp\left(\frac{i 2\pi n t_0}{T}\right)=\sum_{T=1}^{T_{\max}}\cases{1&$T|t_0$\cr 0 &otherwise\cr}=0 \tag{1}$$
which holds iff $t_0$ is not divisible by any integer $1\le T\le T_{max}$. However $t_0$ is always divisible by $1$ so the sum in $(1)$ must always evaluate to something $\ge1$ for any $T_{max}$. In which case we are forced to say

there are no integers $t_0,\,T_{max}$ such that $(1)$ holds and thus:
$$\mathcal{Z}(T_{max})=0, \,\forall T_{max}\in \mathbb{N}.\tag{2}$$

Alternatively, if we chose to exclude the trivial case in which $1$ divides all integers by beginning the outer sum in $(1)$ at $T=2$ we may find a less trivial outcome. Define $t_0$ with $2\le T_{max}\le  t_0\le T_{max}^2$ to be a solution to the equation:
$$\sum_{T=2}^{T_{\max}}\frac{1}{T}\sum_{n=1}^{T}\exp\left(\frac{i 2\pi n t_0}{T}\right)=0$$
then $\mathcal{Z}$ counts the number of integers $t_0 \in \left[T_{max},T_{max}^2\right]$ such that $T\nmid t_0 \forall\, 2\le T \le T_{max}$, in which case:
$$\mathcal{Z}=\pi(T_{max}^2)-\pi(T_{max}) \tag{3}$$
where $\pi(x)$ is the prime counting function.

Proof

If: $$\exists t_0 \,\, \text{such that} \,\,\,T\nmid t_0 \,\forall\, 2\le T \le T_{max},$$ but $t_0$ is not prime, then: $$\exists k>T_{max}\,\, \text{such that} \,\,k|t_0,$$ and:
 $$\exists n=\dfrac{t_0}{k}>T_{max},n\in \mathbb{N},$$
which implies:
$$t_0=nk>T_{max}^2,$$
from which we conclude that $t_0 \in \left[T_{max},T_{max}^2\right]$ is not a solution when $t_0$ is not prime. It is clear that if $t_0$ $\left(2\le T_{max}\le  t_0\le T_{max}^2\right)$ is prime, it is not divisible by $2\le T \le T_{max}$ and is therefore a solution, unless $T_{max}$ is itself prime. From this we conclude that all $t_0$ are primes such that $T_{max}<t_0\le T_{max}^2$ and $(3)$ follows.


MATLAB: Prepare an RGB image for digit and letter recognition

First, many OCR programs automatically perform image processing, including possible color conversion, image filtration, dilation/erosion, and monochrome conversion, in order to improve the quality of OCR operation. Especially good in this are neural-network enabled cloud-based OCR applications.
If you absolutely need local MATLAB program to perform the image processing supposedly preparing the document for recognition, first check if MATLAB or its addons have some kind of image thresholding function: the function which outputs a monochrome (black/white) image given a scanned document. If there is no ready thresholding function, you can try and code it. An introduction into a computer graphics required to successfully code an image preparation operations "for digit and letter recognition" you will find in https://tesseract-ocr.github.io/tessdoc/ImproveQuality and elsewhere, if you search words and phrases highlighted in this answer.
You may find useful this ref, which, although lacking MATLAB code, contains code snippets written in common computer languages.


Applying estimated transfer function to a signal

Your values for m and n are way too big. Matlab starts getting numerical problems for those values above 30, sometimes.

Try reducing them to more reasonable numbers (e.g. 10) and slowly increase them until invfreqz generates bogus outputs.


what is the relation between bandwidth and channel taps

Just as you sample a signal to process or model it digitally, you also sample the channel's impulse response. Each sample corresponds to one tap. When the signal's bandwidth increases, you increase the sampling rate, so your channel model now has more taps.

Usually, a channel is modelled as having a single tap (flat fading) when the channel's impulse response is different than zero only for the first sample. You can then ignore all other taps, since they're zero. This means that the largest path delay through the channel is shorter than the sampling interval. The implication is that a symbol will not interfere with its neighbors.

Note that, if you oversample, you can get a channel impulse response with more than one tap and still have flat fading. Flat fading means that, at the sample corresponding to the symbol time (and all subsequent samples), the channel's impluse response is zero.


Analytically determine a PSD from a transient function

Energy signals, i.e., signals with finite energy

$$\int_{-\infty}^{\infty}|f(t)|^2dt<\infty$$

have zero power and, consequently, a power spectrum that is equal to zero. They do have an energy density spectrum, which is the squared magnitude of their Fourier transform.

You will only get a non-zero power spectrum, according to the definition in your question, for power signals, which have a finite non-zero power defined by

$$\overline{f^2(t)}=\lim_{T\to\infty}\frac{1}{T}\int_{-T/2}^{T/2}|f(t)|^2dt$$

I.e., power signals satisfy

$$0<\overline{f^2(t)}<\infty$$


Can convolution of one signal with different signals give the same answer?

Convolution is filtering.

Consider creating a filter with a rectangular response in the frequency domain.  Then any signal with the same content within the filter's passband, no matter what the signal content outside the passband of the filter, should result in only the content within the passband after convolution with the impulse response of the filter.  Thus, multiple inputs (signals with different stuff outside of the passband) to a convolution with this filter should result in the same output.

(Your opening statement merely swaps the frequency and time domains for the rectangular function, making this relationship harder to recognize.)

This should work even for a non-rectangular and finite length FIR filter with at least one zero in it's transform (maybe at a nice integer submultiple of the sample rate).  Add any magnitude of sinusoidal input at the frequency of the zero, and the convolution result or filter output remains unchanged (except perhaps for numerical issues/clipping/quantization/etc.).


How to prevent distortion after down-sampling?

Read up on the sampling theorem, e.g. https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem

If your feature vectors are still in the time domain and linearly derived from the original audio wave forms, you just need to low pass filter before down sampling.

You need to low-pass filter to avoid aliasing. There many different types of lowpass filters and the best choice depends on the specific requirements of your application (causality, latency, phase distortion, transients, size of transition band, MIPS & Memory, required aliasing rejection, pass band ripple tolerance, etc.)

As you lowpass filter, you lose information. If you need to down sample to a 500 Hz sampling rate all information above 250 Hz (in practice more like 220 Hz) is lost.

Lowpass filters are can be tricky to design depending on your requirements but the code to run them is very straight forward and can easily be done from scratch. Just make sure you use cascaded second order sections and Direct From I or Transposed Form  II.


When calculating the Nyquist frequency should carrier frequency be included

You are correct: a 400 MHz signal needs to be sampled at least at 800 Msamples/s, and that is either impossible or extremely expensive. However, this is not the entire story:

Under some conditions, you can do bandpass sampling, where the sampling frequency is lower than Nyquist. See more at https://en.wikipedia.org/wiki/Undersampling.
You can do analog bandpass filtering and down-conversion to some intermediate frequency (IF), which can then be sampled at reasonable cost/complexity. See https://en.wikipedia.org/wiki/Intermediate_frequency.

Most amateur digital radio receivers I'm aware of use the IF approach, and some also do bandpass sampling. For example, the RTL2832 can only sample at 2.6 Msamples/s, but it can receive signals at up to 2.2 GHz or so. If you see the device's datasheet, you'll see that it supports several different IFs.


First link is broken, the second link is about unrolling loops, so it's not entirely clear (to me) what exactly you mean by "array indexing". Anyway most of this heavily depends on the instruction set of the processor. For example and FIR filter would have to do something like this (in the innner loop)

read filter coefficient
read delayed input sample (state variable)
multiply the two
add to accumulator
decrement coefficient pointer (or decrement coefficient read index)
increment state pointer (or increment state read index)
warp state pointer if needed (these tend to be circular buffers)
decrement loop counter, test for 0
jump to beginning of loop or past the loop depending on previous test

So that's a lot of stuff to do. Especially the "wrapping" part can be quite expensive as it includes a test for wrap around and a conditional which can break the pipeline on many processors.

Anyway on some processors, like a risc, or i86 this can easily take 10-20 cycles, on a dedicated DSP like an Analog devices Sharc or a TI TMS320 these take about half a cycle (can do two in a single CPU clock), so there is a huge dependency on instruction set efficiency. The tricks the second long talks about can help with avoiding the unwarpping part and also reduces the overhead for looping since you only do it once every 4 coefficients.


Stochastic sampling doesn't have anything to do with sampling stochastic waveforms. It simply means that instead of sampling at regular intervals, the waveform is sampled randomly.

Recall that in a sampling scheme per the Nyquist-Shannon sampling theorem, a continuous signal $x(t)$ on $\mathbb{R}$ is sampled as $x[n]=x(nT),\ n\in\mathbb{Z}$, where $T$ is the sampling interval and $f_s=1/T$ is the sampling frequency. If the maximum frequency in the signal is $f_{max}$, then $f_s$ must be such that $f_s\geq 2f_{max}$ so as to avoid aliasing. For ease of comparison with stochastic sampling later on in the answer, let me redefine the sampling in a slightly different form than usual as

$$
\begin{align}
s(t)&=\sum_{n=0}^{f_s\tau -1}\delta(t-nT)\\
x[n]&=x(t)\cdot s(t)
\end{align}
$$
where $\delta(t)$ is the Dirac delta function and $x(t)$ is only sampled on the interval $[0,\tau]$.

If you actually think about it, regular sampling is pretty limiting in practice. Aliasing crops up in several places, and probably a well known and visible effect is the Moiré patterns which can be reproduced at home by taking a photo of regular patterns displayed on a television (examples below).

However, this is always a problem with cameras, but never with your eyes if you were to see the pattern directly! The reason is because the photoreceptors in your retina are not laid out in a regular pattern unlike the CCD in a camera. The idea behind (not necessarily the idea that led to its development) stochastic sampling is very similar to the non-regular layout of photoreceptors in the eye. It is an anti-aliasing technique which works by breaking up the regularity in the sampling.

In stochastic sampling, every point in the signal has a non-zero probability of being sampled (unlike regular sampling where certain sections will never be sampled). A simple uniform stochastic sampling scheme can be implemented over the same interval $[0,\tau]$ as

$$
\begin{align}
s(t)&=\sum_{n=0}^{f_s \tau -1}\delta(t-t_n),\quad t_n\sim \mathcal{U}(0,\tau)\\
x[n]&=x(t)\cdot s(t)
\end{align}
$$

where $\mathcal{U}(0,\tau)$ is the uniform distribution on the interval $[0,\tau]$.

By sampling stochastically, there is no "Nyquist frequency" to talk about, so aliasing will no longer be a problem as before. However, this comes at a price. What you gain in anti-aliasing, you lose by noise in the system. The stochastic sampling introduces high-frequency noise, although for several applications (especially in imaging), aliasing is a much stronger nuisance than noise (e.g., you can see the Moiré patterns easily in the above images, but to a lesser extent the speckle noise).

As far as I know, stochastic sampling schemes are almost always used in spatial sampling (in image processing, computer graphics, array processing, etc.) and sampling in the time domain is still predominantly regular (I'm not sure if people even bother with stochastic sampling in the time domain). There are several different stochastic sampling schemes such as Poisson sampling, jittered sampling, etc., which you can look up if you're interested. For a general, low key introduction to the topic, see

  M. A. Z. Dippé and E. H. Wold, "Antialiasing Through Stochastic Sampling", SIGGRAPH, Vol. 19, No. 5, pp. 69-78, 1985.


$\mathcal Z$-transform if the output is given


There are dozens of publications dealing with shadow-detection, generating shadow-masks, and indeed some that actually remove shadows - such as the ones mentioned in the previous posts. I could add some to the list, if needed. The problem, however, is IMHO far from being solved. For a quick-start, given a shadow-mask, I suggest (and have tried in the past) the two following approaches. They definitely reduce shadows - just not always seamless, and I'm sure there are publications (not by me) dealing with shadow-removal in similar ways.

Gradient Domain Manipulation Techniques as described here (C and Matlab Codes provided):
http://www.umiacs.umd.edu/~aagrawal/ICCV2007Course/index.html
The gradient-integration approach can be used for a number of image-processing problems, see the slides/presentation for further examples.

The general idea:

Compute the spatial derivatives (gradient-images) for all color-channels.
Use the shadow-boundaries from the shadow-mask to generate a weight-mask that is close to zero on shadow-boundaries and increases to one within a specified neighborhood along the shadow-edge, i.e. orthogonal to a given edge-point.
Multiply the weight-mask from (2.) with all the gradient-images to reduce/damp the shadow-boundaries/edges.
Integrate the gradient-images using the code from the above link.
For RGB-images, from my experience, compute the average of the separate channels of the original images, and scale the integrated images to match these values to prevent "funny" color-artifacts.

Brightness manipulation in the original image-domain.

Use the shadow-mask to generate a weight-mask that is one outside shadow regions, has a smooth transition (upwards) across the shadow-boundary and a larger than one scale-factor inside shadow-regions. As suggested in a previous post, the scale-factor can be estimated from the area immediately surrounding a shadow-area using the average-brightness together with the average-brightness of the shadow-region.
Multiply the original images (per channel) with the weight-mask, subject to clipping.

I have also tried using different color-models, e.g. HSV, that present luminance or brightness directly, which can then be modified independent of color (hue/saturation). This works essentially like the brightness-manipulation, i.e. generate a smooth weight-mask and multiply it with the luminance-channel. Maybe the two approaches, i.e. gradient-integration and brightness-manipulation, can be combined in a clever way, but someone has probably tried that before too.

Hope that helps, kind regards, Derik.


How to determine degree of polarization in time domain!

You should probably use one of the algorithms for primary component (PCA) or independent component (ICA) investigation. They heavily relay on covatiance matrices and their eigenvalue-eigenvector decomposition.

If your time series represent an oscillation with dominant direction, the eigenvector for the maximum magnitude eigenvalue would represent that direction. If the signal is some kind of elliptic, the basis of eigenvectors would represent the oscillation ellipsoid dimensions. Degree of polarization is then easily yielded from relations of main-to-other eigenvectors magnitudes.


The transform of a single rectangular pulse is a Sinc function.  Each added translation of that pulse adds a twisted repitition of that Sinc.  The twist will equal 2 pi at all spectral distances that are odd integer multiples of the reciprocal of the repeat period.  Repeat, and eventually all the twisted portions will all cancel out (form circles in the complex plane around zero) except at those points representing the odd harmonics.  Thus an infinite repeating square wave will be left with only non-zero spectrum at odd harmonics in its transform.


I am developing the same algorithm, but I am using C++ one of the biggest  problem I got was that I  did the math (subtraction and sum) using uchar values obtaining black images. Try transforming uchar to float to do the math and then go back to uchar to visualize the image.


The idea is that you can transmit two signals, $I(t)$ and $Q(t)$, over the same bandwidth and at the same time, and still recover each independently of the other. The math is pretty simple. If the transmitted signal is $s(t)=I(t)\cos(2\pi f_ct)-Q(t)\sin(2\pi f_ct)$, then (ignorning noise, ignoring a factor of 1/2 and assuming coherent reception) the receiver can calculate
\begin{align*}
r(t)\cos(2\pi f_ct) &= I(t)\cos^2(2\pi f_ct)-Q(t)\sin(2\pi f_ct)\cos(2\pi f_ct)\\
&=I(t)+I(t)\cos(2\pi 2f_ct)-Q(t)\sin(2\pi 2f_ct).
\end{align*}
After a low-pass filter, the receiver has recovered $I(t)$. Likewise,
\begin{align*}
r(t)\sin(2\pi f_ct) &= I(t)\cos(2\pi f_ct)\sin(2\pi f_ct)-Q(t)\sin^2(2\pi f_ct)\\
&=I(t)\sin(2\pi 2f_ct)-Q(t)+Q(t)\sin(2\pi 2f_ct).
\end{align*}
from which $Q(t)$ can be recovered by an low-pass filtering and multiplication by -1.

The identities used here are: $\cos(A)\cos(B)=0.5(\cos(A+B)+\cos(A-B))$, $\cos(A)\sin(B)=0.5(\sin(A+B)-\sin(A-B))$ and $\sin(A)\sin(B)=0.5(\cos(A-B)-\cos(A+B))$.

Of course, you can also define the transmitted signal as $v(t)=I(t)\cos(2\pi f_ct)+Q(t)\sin(2\pi f_ct)$ and avoid the sign inversion when estimating $Q(t)$. The main reason not to, though, is that the complex envelope of $s(t)$ is $I(t)+jQ(t)$. Since most of the time it is preferrable to work with the complex envelope, $s(t)$ is chosen to produce the most convenient expression for it.


Whenever talking about Cross Correlation of signals it is important to well define it.

If you mean Cross Correlation in most common sense (Signal Processing):

$$ \left( x \star y \right) \left[ n \right] = \sum_{i = 0}^{N - m -1 } {x} \left[ m \right] y \left[ n + m \right] $$

Then of course not.
You can have cases where the DFT of the Cross Correlation has higher absolute values then the DFT of the Auto Correlation.

  The Schwarz inequality answers your question, Yes.
  But don't use pwelch for the cross spectrum. The cross spectrum is complex valued.

The above by @Stanley Pawlukiewicz is wrong.
Simple way to see so is by defining $ y \left[ n \right] = 10 x \left[ n \right] $ which of course mean the cross correlation between those is higher than the cross correlation between $ x $ and itself.

The reason is that Cauchy Schwarz Inequality is about Inner Product Space.
So if we talking about signals as vectors then we need to use the Inner Product which is the cross correlation at $ 0 $.
If we talk about Random Process than we need to use their definition (Which normalize the energy and then everything works).

If you use Cross Spectrum while sticking to the Random Process definition and normalize the Cross Correlation by the variance of the processes than indeed by Cauchy Schwarz Inequality you're guaranteed to see what you got.


Is the modulation on my WiFi different to the modulation on cellular?

From google:
3G-5G appears to use OFDM —along with QPSK og QAM depending on channel conditions.
WiFi 802.11b used some CDMA-dialect, 802.11a and g uses OFDM along with BPSK, QPSK or QAM.
WiFi mainly use frequencies around 2.4Ghz, while cellphones have used lower frequencies around 900Mhz (due to larger distances). Both have recently added capability for using higher frequencies.


Well, yes, everything that you can calculate needs some time to calculate. Viterbi is quite well-behaved with respect to that, since you can decode up to the last couple of bits (only the last constraint length bits are relevant to the current decision), and you can already prepare the path metrics while receiving sequentially, and only need to do the backtracking after reception of the last code bit.
However, in most practical systems, you'd want some sort of interleaving. Also, you'll have filtering and all sorts of operations with inherent delay before you even decode. So, there's group delay and processing latency in your system anyway.
That, together with the fact that Viterbi is not an overly complex method and typically not used with convolutional codes with very large constraint lengths, and that 1 ms is a rather long time for a modern processing system:
You are probably barking up the wrong tree. You haven't defined acceptable latency yet, nor the length of your burst, nor the degree of error correction you need.
Start with that, and figure out (we'll be glad to help) whether processing delays will be a problem or not, before worrying about your decoder.


I figured it out. It turned out that I don't really need that 'Phase Accumulator', and I should reverse the control signal in the NCO. So, I use -PhErr(i-1) as a control signal to NCO. Here is the modified code:

% Siraj Muhammad
% 25/3/2015
% BPSK Demodulator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load RRC.mat

fc = 0.05000001;
phase_offset = pi/7;

N = length(rf_signal);
bb = zeros(1, N);
bb_f = zeros(1, N);
I_f = zeros(1, N);
Q_f = zeros(1, N);
I_r = zeros(1, N);
Q_r = zeros(1, N);
error = zeros(1, N);
PhErr = zeros(1, N);
error_integral = zeros(1, N);

% Loop Filter Coefficients
BW = 100;  % Hz
loop_theta = 2*pi*BW;
C1 = 4*(loop_theta)^2/(1+sqrt(2)*loop_theta+loop_theta^2);
C2 = 2*sqrt(2)*loop_theta/(1+sqrt(2)*loop_theta+loop_theta^2);

for i = 2:N
    % Downconverting to Baseband
    bb(i) = rf_signal(i).*exp(j*2*pi*fc*i+j*phase_offset).*exp(-j*PhErr(i-1));

    % Filtering
    bb_f = filter(RRC, bb(1:i));
    I_f(i) = real(bb_f(i));
    Q_f(i) = imag(bb_f(i));

    % Error
    error(i) = I_f(i).*Q_f(i);

    % Loop Filter
    error_integral(i) = error(i).*C1 + error_integral(i-1);
    PhErr(i) = error(i).*C2 + error_integral(i);
end

figure; subplot 321; plot(real(bb)); title('I Channel')
subplot 322; plot(imag(bb)); title('Q Channel')
subplot 323; plot(I_f); title('I Channel Filter')
subplot 324; plot(Q_f); title('Q Channel Filter')

subplot 325; plot(error); title('Phase Error');
subplot 326; plot(PhErr); title('Loop Filter');

I added a variable phase_offset to test the loop with different offsets, and I can drift the frequency a little bit and the loop will lock very nicely.

I hope this helps somebody as I couldn't find any MATLAB code for Costas Loop online.


How does a sound card translate amplitude values to sound when playing WAV/PCM/FLAC files?

They represent fluctuations from atmospheric pressure.

Positive values = increases in pressure = speaker cone moving outward
Negative values = decreases in pressure = speaker cone moving inward
Zero = no change in pressure = speaker cone in its relaxed center position = 1 atm

Since we can't hear a constant pressure, the constant pressure is averaged out and ignored, and only the relative vibrations are important.


How to generate AWGN with given power spectral density

Here is a Matlab snippet that might help. You should be able to apply the result to your favorite PSD routine to see the dbm/Hz in the frequency domain.

% create complex noise with a desired dBm/Hz
fs=1e6;             % whatever you want
N=1e5;              % # samples of noise
R=50;               % ohms
dbm_hz = -135;      % goal
bw = fs;            % complex noise with BW = sample rate
dbm = dbm_hz + 10*log10(bw);
dbw = dbm - 30;
watts = 10^(dbw/10);
vrms = sqrt(watts * R);
n = randn(N,1) + j*randn(N,1);
n = vrms * n / std(n);
meas_dbm_hz = 10*log10(sqrt(mean(n.*conj(n)))^2 / R) - 10*log10(bw) + 30


The difference of gaussian (DOG) is the convolution of input image by difference of two gaussians usually with different standard devitations($\sigma$). The basic idea behind this is to capture edges or gradients in the images that are simplified by the gaussian with larger $\sigma$  but preserved by the smaller gaussian. This is form of scale space derivative, which basically analyses different frequency components in the image, by changing the $\sigma$s  of gaussians. Thus, its like a band-pass filter that rejects certain frequencies in the input.

If we consider an constant image with its value centered around the means of the gaussians, the result of the DOG would be a zero - indicating no gradient.Try this out for a basic tutorial.


I am not an Octave user, but the Video package appears to do these things.

Good luck!


In "ordinary" compressive sensing you take measurements $y = Ax$ and try to reconstruct $x$ from $y$ where $A$ - of dimension $(M \times N)$ - has been generated with no knowledge of x, except that it is sparse. $A$ is typically generated randomly.

In adaptive compressed sensing, you in principle take each individual measurement $y_m,\ m = 0, 1,\ldots M$, trying to adapt you next measurement $y_m$ cleverly so you get the best possible measurement for reconstruction of $x$, given what you have already observed in $y_0, y_1,\ldots, y_{m-1}$. Effectively, this means determining the rows $A_m$ of $A$ one at a time, taking into account the measurements that previous rows $A_0, A_1,\ldots, A_{m-1}$ gave you.

One example of adaptive compressive sensing could be Malloy & Nowak, 2013.


One representative time series signal from multiple time series data

Decision of representative is usually formulated as an optimization problem. The optimal representative minimizes sum of some dissimilarity measures from the observed signals. The average you indicated is optimal in the sense of minimization of Euclidean distance. If you want to minimize feet of perpendiculars, you can use the first principal component.

Now the problem is how to choose the dissimilarity measure. In signal processing, we often consider some probabilistic model and use the dissimilarity measure proportional to the negative log likelihood. This is often referred to as maximum likelihood estimation. For example, Euclidean distance corresponds to the negative log likelihood of mean parameter of Gaussian distribution. Even using the same Gaussian distribution, there can be various ways of modeling, and various dissimilarity measures are derived such as Mahalanobis distance and Itakura-Saito distance.