How can a signal be imaginary?

It's not imaginary in the everyday sense of meaning something that doesn't exist.
"Real" and "imaginary" have technical meanings in signal processing, and more widely in mathematics, that are different from their everyday meanings.
In signal processing "real" denotes a signal component that is in phase with (i.e. has a 0° phase shift from) some reference signal.  Imaginary denotes a signal component that is in quadrature with (i.e. has a 90° phase shift from) the same reference signal. The reference signal can come from a local oscillator. (In DSP or SDR equipment the local oscillator might be a mathematical representation of one, rather than an analog circuit.) In a receiver we might seek to synchronise the reference signal to the received signal.
In signal processing equipment we also refer to real and imaginary as I and Q for in-phase and quadrature, respectively. Some equipment has sections which are divided into separate I and Q channels. You can measure the signals flowing through a Q (imaginary) channel, implemented in hardware, and find they are just as real (in the everyday sense) as the signals flowing through the I (real) channel.
The term "imaginary" was originally coined in the 17th century by René Descartes as a derogatory term and has been confusing students ever since.
For why $e^{j \omega t} = \cos  \omega t + j \sin \omega t $ see Euler's formula. Mathematicians use $i$ to represent the imaginary unit. Electrical engineers use $j$ to represent the imaginary unit because we use $i$ to represent (instantaneous) current.  You can use either representation but please choose one and stick with it.
(This question was originally posted in the amateur radio SE group. This answer was drafted as a reply there and leans towards radio usage.)


Does dividing the magnitude spectrum of white noise by sqrt(2) give an RMS magnitude spectrum?

However, an FFT of a noise signal indicates magnitude in bins at equally spaced frequencies, which means the magnitude represented in each bin, is the magnitude of a periodic signal (sinusoid).

No; it just means that when projecting the noise signal onto a sinusoidal signal, you don't get nothing.

Assuming your noise realization in each is independent identically distributed (usually the way we model noise!), then the (rotated) sum of such noise samples approximates a normal distribution, no matter how the (finite-variance) original distribution looked like. The magnitude in each bin is an approximation (because by definition, your DFT has finite length!) to the noise power density; in an uncorrelated noise process, that's the variance==power==rms².


I trigger recording when power of current sample is X dB greater than last Y samples, what is the terminology here?

"Running average" is fine. You can also call it the moving average.

I suppose you could call the variable associated with the 10dB number the signal detection threshold.


How can $a_n$ for the trigonometric Fourier series be $0$ for a non-odd square signal?

Indeed, if the cosine coefficients ($a_n$) of the Fourier series of a function $s(t)$ are all zero, then the function is the sum of sine terms, which are all odd, implying that the function $s(t)$ is odd.

In your case, however, not all $a_n$ are zero; in particular, $a_0=2$. This means that the wave has a "DC" or mean value of $a_0/2=1$. Since $a_n=0$ for $n>0$, then the function $$u(t)=s(t)-a_0/2$$ is indeed odd, and now $a_n=0$ for all $n$.


I am afraid there's not much you can do.  The voice part seems to have gone through the equivalent of a low pass filter with a cut off of around 1000 Hz.  Basically, all of the speech components above 1000Hz are gone.

The filtering action may not have been an intentional filter, but may have been due to the improper biasing of the tape during recording.  If it is an old tape, it may simply have deteriorated over time.  Also, the playback head may have needed degaussing.

Running it through a high pass flattens the frequency response, but doing so pushes the level down so far it drowns the signal in the noise - no help there.

The best result I got was from using a very steep low pass with a cutoff of 1000Hz together with a very steep high pass with a cutoff of 160 Hz.  That gets rid of the noise by only passsing what is left of the actual speech, but obviously it can't recover what was lost.

Your real problem is not the noise, it is the lost frequency range.

This is the spectrum of the bad recording:

This is the spectrum of the good recording:

As you can see, there's a lot of stuff missing from the bad recording.  So, it isn't simply a problem of removing noise.  The problem is that there's stuff that's just GONE.

Look at the range from 1000Hz to 7000Hz.  There's lots of stuff there in the good recording, but in the bad one it is just a flat spectrum a good 30dB below the voice peaks aroudn 400Hz.

Some of what's missing might be buried in the noise, but recovering it would cause artifacts that are worse than the noise and muffled sound.

Looking at just the noise, it doesn't seem like there would be much to recover out of it.  It looks just like the portions with speech (except for between 160Hz to 1000Hz,) so anything that is in there is going to be buried really deep.


As pointed out in the comments, the inverse system is not realizable. However, it is possible to write down an expression for the impulse response of the inverse system. We want a system that results in a Dirac impulse when convolved with the given rectangular pulse $h(t)$. This can be done as follows. We delay the rectangular impulse by $T/2$, so it starts at $t=0$. Then we add an infinite number of delayed versions of it such that the result becomes a step function $u(t)$. Each delay is an integer multiple of $T$:

$$u(t)=\sum_{n=0}^{\infty}h(t-T/2-nT)\tag{1}$$

The system that achieves this has an impulse response given by

$$f(t)=\sum_{n=0}^{\infty}\delta(t-T/2-nT)\tag{2}$$

Now we have a step function, from which we can generate a Dirac impulse by taking the derivative. So the inverse system is given by the concatenation of the system with impulse response $f(t)$ given by $(2)$ and an ideal differentiator. The total impulse response of the inverse system is then

$$g(t)=(f\star\delta')(t)=\sum_{n=0}^{\infty}\delta'(t-T/2-nT)\tag{3}$$

where $\delta'(t)$ is the distributional derivative of the Dirac impulse.

Clearly, the system with impulse response $g(t)$ given by $(3)$ is not realizable, but as a theoretical concept it is interesting, just like other ideal systems (e.g., ideal frequency selective filters, differentiators, Hilbert transformers, etc.).


The step response you have is non-minimum phase. That means various things:

There are zeros (or poles, if the system is unstable) outside the unit circle (in discrete time).
The peak of the impulse response is delayed from the origin.

Below is an example showing an "ideal" delayed step, a non-minimum phase filter impulse response and the output.

As you can see, the output exhibits similar characteristics to what you have plotted.


Is it possible to approximate given skewness and kurtosis values using an IFFT?

When you know that the process that generated the data is non-linear (and you are lucky to be in full control of the acquisition), you can try Attractor Reconstruction.

This technique attempts to reconstruct the trajectory a system may be taking through phase space, that results in the complex signal that is being recorded. The signal itself, may be showing chaotic behaviour with seemingly unpredictable spikes or a trace that looks like broadband noise. However, through Taken's theorem it is possible to find a "delay embedding" within which this behaviour is "explainable".

In practice, this looks like two (for example) pointers over the signal, separated by some distance (or delay) $k$. Given a one dimensional signal $x[n]$, this results in a two dimensional signal composed of $x[n], x[n+k]$ with $n,k \in \mathbb{N}, k>0$. If $k$ (and the number of dimensions, here 2) is chosen correctly, the subsets of points (here $x[n], x[n+k]$) as plotted in a Cartesian plain, would reveal an attractor.

So, given this analyical framework, the task of "synthesizing a signal" becomes that of recovering the attractor, capturing it in a model and then playing it back at will, to generate the data.

As you may suspect, the choice of $k$ and number of dimensions is the main headache here. Theoretically, this should be equal to the fractal dimension of the attractor (which is not always known) but there are many different metrics and methods for estimating that.

An example of the whole process is available in this paper or, much more extensively in Kantz's book.

Finally, in more recent years, that "modeling" part of the attractor from above could not have escaped the Machine Learning treatment, so for the more up to date approach, you could refer to something like this one or even a GAN approach. These are everywhere these days.

You can see the...attraction here, especially when you are in control of the experimental conditions. I think that something like this would capture the dynamics you are after in a better way and manage to reproduce the conditions you are (were (?)) facing.

Hope this helps.


Trying to understand downsampling and then upsampling

Downsampling can be defined as follows:

$s_{out}[n]= s_{in}[n\cdot M]$ (taking every M-th signal sample)

After applying Z-transform we obtain:

$\bar{s}_{out}[z]=\dfrac{1}{M}\sum_{m=0}^{M-1}\bar{s}_{in}\left[ z^{1/M} \cdot w_{M}^{m} \right]$, where $w_{M}=e^{-2\pi i /M}$

By substituting $z=e^{2\pi i f_n}$ ($f_n $ is the normalised frequency: $f_n = \dfrac{f}{f_s} \in [0, 0.5]$

This yields following relationship between input and output spectrum:

$\hat{s}_{out}[f_n]= \dfrac{1}{M}\sum_{m=0}^{M-1}\hat{s}_{in}\left[ (f_n-m)/M \right]$

You can understand that as a fact that your spectrum is being copied M times (even outside of your original sampling frequency) and stretched by factor of M.

Upsampling on the other hand is governed by following equations:

$s_{out}[n]=\begin{cases} s_{in}(n/M) & \text{n} \; mod \;\text{M} =0 \\ 0 & \text{otherwise}\end{cases}$

After applying Z-transform:

$ \sum_{n} s_{out}[n]\cdot z^{-n} = \sum_{n} s_{in}[n]\cdot z^{-nM}$, which gives: $\bar{s}_{out}[z]=\bar{s}_{in}[z^M] $

By doing same substitution of a unit circle instead of z variable:

$\hat{s}_{out}[f_n] = \hat{s}_{in}[M\cdot f_n] $

You can imagine that this equation is squeezing your spectrum by factor of M, bringing at the same time whatever is outside of your sampling frequency range.

Plots you've obtained are great example how it works!

P.S. There is a dft tag.


Empirical Mode Decomposition and Sparsity

The idea of being sparse in a certain basis can be extended and one can talk about being sparse in a frame, or even, being sparse in a dictionary. Now we may lose orthogonality, have redundant signals in our dictionary and no longer have unique signal expansions; which is OK if we only care about sparsity.

The EMD expansion is sparse in the sense that you can express a signal as a linear combination of just a few IMFs.


How to increase the spectral resolution?

You seem to have a good grasp of the tradeoffs. When using short-time Fourier analysis like you are, there is a version of the uncertainty principle at play. Increasing your time resolution (in your case, using a shorter DFT) results in coarser frequency resolution, and vice versa. That is, the time-bandwidth product is a constant.

The way to increase your STFT's spectral resolution is to increase the duration of time that the transform covers, as you noted. If you truly need to be able to resolve frequencies that precisely (within a few Hz of one another), then you need to observe them for a long enough period of time to discern them. If you know a priori some characteristics of your signal, and conditions are favorable (i.e. SNR is high enough), then you might be able to get the job done with coarser resolution (and therefore a shorter transform).

For instance, if you know that your signal is likely to be a single tone somewhere in a particular band, and you want to know its frequency precisely, then you don't necessarily need to use a really long DFT. Instead, you can use a shorter DFT, then use peak interpolation techniques to give a sub-bin estimate of where the peak actually lies.


How many non-zeros samples in convolution

Actually, you can't say how many non-zero values $a[n]$ is going to have because, depending on the values that $b[n]$ and $c[n]$ take, $a[n]$ may equal 0 for some values of $n$ or not.

Nevertheless, I think that you want to know the length of the sequence $a[n]$. As a general rule, for two sequences of length $P$ and $L$, the linear convolution between them returns a sequence of length $N=P+L-1$. In this case, the maximum length of $a[n]$ would be $N=30+40-1=69$ (it can be less than 69 if $b[n]$ and $c[n]$ are such that the convolution between them returns null values for the first or the last values of $a[n]$).

The reason why the resulting sequence has a length of $N=P+L-1$ can be easily seen by doing the convolution graphically. Let's assume for simplicity (although the reasoning without these assumptions would be exactly the same) that the first non-zero value of $b[n]$ appears at $n=0$ and that the same thing happens with $c[n]$. Let's remember what the definition of convolution is:

$$a[n]=b[n]*c[n]=\sum\limits_{m=-\infty}^\infty b[m]c[n-m]$$

This means that we should have in mind $b[m]$ and $c[-m]$, and then start moving $c[-m]$ from $-\infty$ until we get our first non-zero value for $a[n]$. This happens for $n=0$. Then the two sequences always return some information until the last value of $b[m]$ is in the same spot as the first value of $c[n-m]$. This happens for $n=L+P-2=68$. So that leaves us with the length I wrote before: $N=69$.


Is there an intuitive explanation for what complex basebanding a signal is?

What exactly is the complex baseband of a signal? I sort of understand what "baseband" means, but where does the complex phase come from?

A baseband signal is one where the center frequency is $\approx 0$ Hertz. There may be some small offset, but the idea is that it is no where near a passband signal whose center frequency is in the millions or billions of Hertz.
For a real signal, the frequency spectrum is symmetric. Due to the up-conversion process (converting baseband to passband so that the center frequency is no longer $0$ Hertz), we can put data on two orthogonal carriers, $\text{cos}$ and $\text{sin}$, in the passband signal:
\begin{align}
x(t)&=\text{Re}\bigg(s(t)e^{j2\pi f_ct} \bigg) \\
&= x_I(t)\text{cos}\big(2\pi f_ct\big) - x_Q(t)\text{sin}\big(2\pi f_ct \big)
\end{align}
where $s(t)=x_I(t) + jx_Q(t)$, is the complex baseband signal. For a nice explanation of what is probably happening in your MATLAB function on how to convert passband to baseband, check out this previous answer: https://dsp.stackexchange.com/a/43281/31316.

Related, does the phase have any relation to the signal? If I input a sinusoid, what should I expect the phase output to be?

There is a relation...it is the phase, that is all. The phase is helpful though and can help to pack more bits of information into a single symbol.


How to calculate dB of a transmitted power

If you have a system that causes 80% of its input power to be lost, then 20% of the input power will remain at the output. The power gain of the network is therefore:

$$
G = \frac{P_{out}}{P_{in}} =0.2
$$
$$
G|_{dB} = 10 \log_{10}(G) = 10 \log_{10}(0.2) \approx -7 \text{ dB}
$$


Finding frequency respone of a differential/integral LTI system

Frequency response deals with the steady-state behavior. Differentiating with respect to $t$ yields:

$$\frac{d^2y}{dt^2}+2\frac{dy}{dt}+3y(t)=\frac{dx}{dt}-5\frac{d^2x}{dt^2}$$

The Laplace transform is $$(s^2+2s+3)Y(s)=(s-5s^2)X(s)$$
and we have

$$H(s)=\frac{Y(s)}{X(s)}=\frac{s-5s^2}{s^2+2s+3}$$
I think you can easily do the rest.


Detect ending of frame in H.264 video stream

I think Technically speaking the only way to detect when a frame ends is to see where the next one begins. That means, you should search for the next start prefix code 0x000001. If you go to the previous data before this (which in the general case is the byte alignment syntax elements) and add your RS codes, a compliant decoder should ignore all that until the prefix code is seen.

As for the end of GOP, I'm afraid that the GOP concept does not exist in the Video standard, in theory you could start receiving frames which use previous frames as reference, but there is never a repeating structure. You can detect the end of a coded video sequence, but that does not necessarily coincide with what we know in the video community as "the GOP".

Edit:
In reply to the comment:
Lets suppose the bit stream is this (I am going to invent, it will be all garbage):
0x 00 00 01 e5 d8 31 04 ab 80 00 00 01 8e 56 04 94 cd 80 10 00 00 01 ...

You know for example that 00 00 01 e5 d8 31 04 ab 80 is the first frame, so you take all those bits and generate a CRC for them (or any other error correcting technique that you want), lets suppose you insert one byte of CRC ( that hypothetically will allow you to detect up to 5 bits in error and correct 3 for example), then your new bit stream would be

0x 00 00 01 e5 d8 31 04 ab 80 fc 00 00 01 8e 56 04 94 cd 80 53 10 00 00 01

here the CRC would be the bytes in boldface.

As for the GOP. Remember that even though when you decode you can see a repeating structure of frames "it is not required" to be like that forever. IF you want you could count the number of frames between I frames and in the practical cases it will work. Most H.264 encoders use a structure like I  {[B] (xN) P} (xM) I ...

Examples:

I B B P B B P B B I B B P B B P B B I B B P B B P B B ... etc

I B B B P B B B P B B B I B B B P B B B P B B B I B B B P B B B P B B B ... etc

But if you start counting and see in first example that the number of frames between two I frames is 8, nothing guarantees that in the next cycle it would be the same. Usually it will be, but there is nothing in the standard saying that you have to do it like that.

Hope it helps


How to compare 2 matrices of MFCCs, keeping the time information

You are misunderstanding how DTW works. You can actually adjust the frame deletion/insertion costs to adjust how much temporal mismatch you can tolerate.

We would have: DTW("Puuuum pum tchak", "Puuum pum tchak") < DTW("Puuuum pum tchak", "Pum pum tchk") < DTW("Puuuum pum tchak", "Purr purr purr").

That is to say, if there is a sample with a perfectly good timing, it'll always yield a higher DTW score than a sample with the wrong timing. The trade-off between similar timbre and similar timing can be adjusted by tweaking the DTW alignment costs. One useful thing with DTW is that it finds the optimal alignment between the "query" and the found sample - so you can use the path it founds to timestretch/phase-vocode the result and align it perfectly with the query.

It is not clear what kind of sensitivity on timing your application has; which features matter the most (overall timbre? onomatopoeia which is pronounced?); and what kind of processing you could do to the found clips to align them even better with the query. If rhythm/timing is absolutely essential, you might want to process each sample by an onset detector ; and restrict your search to the set of clips which have the right number of onsets, and the correct inter-onset times.

If you have absolutely no tolerance on differences in rhythm, then you could very well compute the average of pairwise distances between the two clips (padded with zeros if their length does not match). That would be a bit silly, but your application might want that.

Finally, if both your query and your samples can accept the same domain-specific representation (say a beatboxed rhythm and a drum loop can both be transcribed into "bass drum" and "snare drum" events) - it might help to transcribe all clips and to match the transcriptions instead of the audio signals.


This answer is written in terms of low-pass but  the ideas directly translate to high pass.

A butterworth filter (low pass) of order $N$ with cutoff $\Omega_c$ is one which has a maximally flat passband, and has the squared magnitude response $|H(j \Omega)|^2 = \frac{1}{1+\left(\frac{j \Omega}{j \Omega_{cutoff}}\right)^{2N}}$. [I'm just writing it this way so that if you try to get a realizable transfer function via spectral factorization it is easy].

So, you specify the cutoff $\Omega_c$ and it gets sharper as you increase the filter order. Strictly speaking, from the the squared magnitude response, you see that $|H(j \Omega)|^2=1$ if and only if $\Omega=0$, i.e. only the DC component is un-modified. However, from that formula, you can calculate things like frequencies where the magnitude is within say 99% of the original magnitude, by finding $\Omega$ such that $|H(j \Omega)|^2 \geq 0.99^2$ (note that $|H(j \Omega)|^2 $ is monotone decreasing in $\Omega$, so your range will be of the form $[-\Omega_{0.99},\Omega_{0.99}]$). Similar properties of other filters will hold as well - see Appendix B of Discrete-Time Signal Processing 2e by Oppenheim, Schafer and Buck for more details. Make sure you're familiar with the basic types of filters as well -- elliptical, butterworth and chebyshev to make sure your characteristics you desire are matched to what you want.

If you want to design filters with particular passband and stopband characteristics in the FIR case (such as controlling the maximum error), you may want to look at the Parks-Mclellan algorithm (described in Discrete-Time Signal Processing 2e by Oppenheim Schafer and Buck, ch. 7).

As for Matlab usage, you can plot the magnitude response of the filter and determine these things. It can be done from the filter visualization tool. TI also makes a filter design tool which is quite nice, called FilterPro.


What is the purpose of dividing a signal into subbands when compressing an audio file?

We don't hear all frequencies equally well. Therefore, if we split the frequencies into subbands, we can give more bits to bands which we hear well, and less bits to bands we don't hear that well (which is determined by things like doing listening tests). Its usually easier easier to control the error of the compressed signal on a subband level than on the whole signal level.

This will typically give us less distortion for a given data rate than trying to encode all frequencies with equal number of bits.

An analogous experiment you can do is with (natural) images: Take the DFT or DCT of an image, and throw out the higher frequencies (i.e. high in the x,y directions). You still get a pretty decent approximation if you keep only a few low frequencies. So, instead of throwing them out, imagine you quantized them -- you'd want to give less bits to the higher frequencies than the lower frequencies, since you're not getting much perceptual gain by throwing lots of bits in the higher frequencies, but you are getting more in the lower frequencies. (This underlies JPEG; See something like Gonzalez and Woods' Digital Image Processing for details on how exactly JPEG works).


How to compute correct phase in FFT even after applying phase unwrap and zeroing round off error?

There's an error in your script. Your sampling frequency is not 100 Hz like you intended but is actually 99.9 Hz.  Therefore you don't actually have a whole number of cycles in your FFT which creates some leakage. This leakage will also impact the phase estimation.
If you want 1000 points at 100 Hz, the last point should not end at 10 seconds but at 10 seconds - 1/100 = 9.99. You should be careful when using linspace to define a time vector.
Use this n = 1 : 1000; t = (n-1)/100;


The answer is no:

See the following rearrangements of numbers from 1 to 81.
In the left  case the 3x3 median of 3x3 medians will be

median([77, 72, 67, 62, 57, 52, 47, 37]) = 57,

in right case
median

([5, 15, 23, 32, 41, 50, 59, 68, 77]) = 41

While the median of the whole 9x9 block is 41 for both cases.


Designing A Filter Using SPTool In Matlab

John Smith,

You are designing a low pass filter, which means that your filter will be applied, (convolved) against a signal, whereby it will attenuate certain frequencies, and leave other frequencies more or less untouched.

So in a low pass filter, generally that means that you will filter out (heavily attenuate) frequencies that are 'high'. Your $f_{pass}$ just quantifies this for you. How high is high? Well, if you want a 'cutoff' at 200 Hz, then that generally means you want to keep frequencies at roughly 200 Hz and below, so pick $f_{pass}$ to be 200 Hz.

But bear in mind, you usually have a specification about transition bandwidth. If your filter is infinitely sharp, you can say that you accept all frequencies less than 200, and attenuate all frequencies above 200. Realistically however you might want to accept all frequencies less than say, 130, and then from 130 to 200 you are transitioning, and then heavily attentuate all frequencies greater than 200.

Regarding the sampling rate, that is really dependent on your simulation. If you pick 10,000 Hz as your sampling rate, then you should make sure that when you use your filter on a signal, that it too was sampled at 10,000 Hz.

It is good that you are using MATLAB. It is a very powerful tool in understanding DSP concepts.


You can only "ignore Nyquist" if you a-priori know that there is no higher frequency content (above your desired or required noise floor) that will alias down (or up!) to your band of interest.  But since your signal has content at and above half your sample rate, you will need to low pass filter that signal before sampling.

You will also need a sampling rate ABOVE 1000 Hz, if you are interested in a finite length of signal at 500 Hz.


Output of Zero-forcing equalizer doesn't match theoretical calculations

It seems like there might be something weird going on with that blog, but I only quickly looked. Here is the idea for the ZF equalizer:

Figure out the convolution matrix.

convMatrix = convmtx(h, length(x));

This will be a $9 \times 7$ matrix because the input vector $\mathbf{x}$ is length 7 and the output convolution will be of length $7+3-1=9$. It is not a square matrix but we can still get an inverse by computing the Moore-Penrose inverse (https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse), which in MATLAB can be done by calling:

invConvMatrix = pinv(convMatrix);.

Choose your delay. invConvMatrix is a $7 \times 9$ matrix and we are multiplying with an all zero vector except for a single one at some delay $\mathbf{\delta}$. Lets choose the delay to be 5, so we have:
$\mathbf{\delta}=[0, 0, 0, 0, 1, 0, 0, 0, 0]^T$.
Compute the ZF equalizer:

 zfEqualizer = invConvMatrix * d;

Generate the received signal, and generate signal after equalization

 y = conv(h, x);
  z = conv(zfEqualizer, y);

Now lets take a look at the resulting signals, and specifically notice where the equalized signal starts, right at the delay of 5 that we choose earlier. I hope this helps.


Visible light communication and LED modulation

OFDM in general is a complex baseband technology, which means it produces an I and a Q signal, which have to be mixed with a carrier wave and it's 90° shifted version, and added up before it can be used to excite a real-valued emitter (like an antenna or an LED driver).

By "I made the IFFT output real" you probably mean that you supplied the IFFT with a symmetric signal? That's nice, because it circumvents the need for a carrier frequency, but it also inherently halves your available bandwidth! Also, I can't fight the feeling that the discrete Fourier transform is not really what you're looking for if you just want to produce a real signal – the discrete cosine transform might be better suited for you.

The core question here is why you want to use OFDM in the first place.

The main motivation for OFDM comes from something that is highly inherent to the RF channel: Channels of high bandwidth aren't flat, so the equalizer for something like 20MHz of WiFi channel gets really ugly. Here, you can gain a lot by dividing your channel into orthogonal subchannels, and equalize them with a much simpler equalizer individually. I wasn't aware the same problem applies to the monochromatic VLC channel!


Custom ringtone sounds good on computer, can't be heard on phone

Plug your phone into some large speakers.  If you can hear the tone much better, then the problem might be the small size of the speakers on the iPhone, which roll off low frequencies.

To prevent this, make sure any low pitched notes in your ringtone are made with timbres rich in strong higher overtones.

Added: If you create a ringtone from a stereo file, it's probably best if both channels contain exactly the same data.


Filter design for an unusual EEG experiment

Is my filter "reasonable"?

I would say no. The time domain ringing time of your filter far exceeds the length of your signal, there isn't enough signal for the output to stabilize.

Can I do better?

Difficult to say. Filter design involves complicated tradeoffs between many different factors and the "best" really depends on the specific needs & sensitivities of your application. You have partially explained what you want to do but not in enough details to give a recommendation. I'm skeptical about "looking at the phase" approach. Signal phase has only meaning if you have a well defined time reference for $t=0$, which most signals don't.

Can you recommend to me some learning materials for a filter design.

I personally like that one https://ccrma.stanford.edu/~jos/filters/ but I'm also mostly and audio guy.

(*) This is obviously cut from continuous EEG signal - can it be just cut (boxcar window),
or should I apply something different (eg. Hamming window)?

Depends a bit on your application details. Ideally your filter "ringing" is significantly shorter than your signal. Then you can just take a "large" chunk of signal, filter it and cut off the edges to remove the filter transients.

(*) Is my expectation justified?

Can't say without looking at your signal first. A sufficiently narrow band signal will look like an amplitude and/or phase modulated sine wave. However, if your original signal is very noise-like, you'd have to make it VERY narrow to look like that.


How to derive stability of analogue and digital systems

In the discrete-time case you can use the $\mathcal{Z}$-transform:

$$\mathcal{Z}\{x(n)\}=X(z)=\sum_{n=0}^{\infty}x(n)z^{-n}\tag{1}$$

Equation (1) is the unilateral $\mathcal{Z}$-transform which is useful for causal signals which are zero for $n<0$. Now consider the signal $x(n)=a^n$, $n\ge 0$:

$$X(z)=\sum_{n=0}^{\infty}a^nz^{-n}$$

This series converges for $|a|<1$ and in this case $X(z)$ is given by:

$$X(z)=\frac{1}{1-az^{-1}}=\frac{z}{z-a}\tag{2}$$

From (2) you can see that $a$ is the pole of $X(z)$, and for $x(n)$ to be stable (if seen as a system's impulse response), and for $X(z)$ to exist, $|a|<1$ must hold, i.e. the pole must be inside the unit circle of the complex $z$-plane. This example is analogous to the continuous-time example in your question.

Also take a look at the concept of BIBO-stability (bounded-input-bounded-output stability), which leads to the same conclusion: for continuous time signals/systems the poles must be in the strict left half of the $s$-plane, and for discrete-time signals/systems the poles must be strictly inside the unit circle of the $z$-plane. The condition on the time-domain signals is that they be absolutely integrable in the continuous case, and absolutely summable in the discrete-time case:

$$\int_{-\infty}^{\infty}|x(t)|dt <\infty\\
\sum_{n=-\infty}^{\infty}|x(n)| <\infty$$


How to calculate the energy levels of frequency intervals of an audio signal?

You may research for spectrogram. The above comments from pichennets and B Z are correct, but you have to pay attention to the effect of just take some subset of coefficients and calculate the power. When you take this way you're applying a retangular window in frequency, and some peaks in the frequency domain power. You must apply a window in time domain before take the FFT and calculate the power.

I want to add this as a comment, since it's not a real answer but I'm not abble to do this.


BER Analysis of Dual hop Decode and Forward Relay

I am working with a similar kind of network. I don't know if it is correct, but I think the most proper way is to consider the SNR from the first transmitter, not considering what the relay do to the signal. This because you want to analyze the performance of the end-to-end communication. SO you send the signal with some energy, and you want to know the performance of your signal once you have sent it with that energy.

That's the way I am doing, and comparing it to the direct transmission, I am getting a very similar performance: a difference of less than 1dB for the same BER.

I am not amplifying the signal in the relay, so I am just doing the decode and forward (DF). I think it is reasonable to get this "worst" performance in the DF relay compared to the direct transmission because you have two channels, therefore you have two components of noise. You perform the DF in the relay, but it can decode wrong.

Let me know if this answer was useful, and if you have any more questions. I am not an expert in this area, but I am working on it.


frequency characterritic of activ high pass filter

There are many online calculators with plot etc. options you could use. An example: http://sim.okawa-denshi.jp/en/Fkeisan.htm


Neural network with much less positive examples

There are a couple of things you can do to treat the class imbalance.

A way is to use an SVM-loss (Hinge-loss) function in the network, where the margin around the decision boundary is maximized. The good part of this: Because this function considers the support vectors only, the data points which do not lie close to the margin will not be taken into account. Of course, with more data, you almost guarantee to change the decision boundary and therefore the support vectors. However, especially if your data is spread around the boundary (if you can sample some positive-like true negatives), then you should get close to the correct minimum.
This is an easy way to get around without much effort.
You could augment (not just duplicate) the positives more than you augment the negatives, so that the data becomes balanced.
The weighting of the samples is already discussed, and I think it is valid.


Calculate Phase angle between two signals i.e. digital Phase meter

A quick and dirty method is to compare the energy of the sum (L+R) to the energy of the difference (L-R). This gives a a crude single-number phase difference and is good at finding things like polarity flips etc.

As Dilip pointed out, the concept of phase difference between uncorrelated signals is meaningless. Even if they are correlated, the phase difference tends to be a function of frequency, so a single number isn't particularly helpful either. Most real world stereo audio signals are partially correlated, which complicates the matter further. In order to really answer the question, you need to properly define what you mean by phase difference and what you are planning to do with it.


Why do we need to have unity gain at DC for a lowpass filter?

There is no mathematical explanation and it's typically just a convenience for most applications. If you scale so that it's unity at 0 Hz, unity is also the maximum gain at any other frequency so you are less likely to clip or overdrive something in the signal chain (it's still possible though).

It's often based on the assumption or requirement that "frequencies below the cutoff should be unaffected" by a low pass filter. The closest way to get this is to have unity gain at DC.


What is the peak resonance of convolving with a sine FIR filter?

If your signal was a windowed complex exponential, then the maximum of the magnitude of its discrete-time Fourier transform (DTFT) would equal the frequency of the complex exponential:

$$x[n]=\begin{cases}e^{jn\omega_0},\quad&n\in[0,N-1]\\0,\quad& \textrm{otherwise}\end{cases}\tag{1}$$

The DTFT of $(1)$ is easily computed from the convolution of the DTFT of a rectangular window $w[n]=1$ for $n\in[0,N-1]$ and zero otherwise, and the DTFT of $e^{jn\omega_0}$:

$$\textrm{DTFT}\{w[n]\}=W(e^{j\omega})=e^{-j(N-1)\omega/2}\frac{\sin\left(N\omega/2\right)}{\sin(\omega/2)}\tag{2}$$

$$\textrm{DTFT}\{e^{jn\omega_0}\}=2\pi\delta(\omega-\omega_0),\qquad |\omega|<\pi,\quad |\omega_0|<\pi\tag{3}$$

Note that the DTFT in $(3)$ is of course $2\pi$-periodic.

With $(2)$ and $(3)$ we can write the DTFT of $(1)$ as follows:

$$\begin{align}X(e^{j\omega})&=W(e^{j\omega})\star\delta(\omega-\omega_0)\\&=W\left(e^{j(\omega-\omega_0)}\right)\\&=e^{-j(N-1)(\omega-\omega_0)/2}\frac{\sin\left(N(\omega-\omega_0)/2\right)}{\sin((\omega-\omega_0)/2)}\tag{4}\end{align}$$

The magnitude of $(4)$ clearly has a maximum at $\omega=\omega_0$ $(+2\pi k)$

However, for a windowed sinusoidal signal we get two contributions of the form $(2)$ centered at $\pm\omega_0$ because $\sin(n\omega_0)=(e^{jn\omega_0}-e^{-jn\omega_0})/2j$:

$$y[n]=\begin{cases}\sin(n\omega_0),\quad&n\in[0,N-1]\\0,\quad& \textrm{otherwise}\end{cases}\tag{5}$$

$$Y(e^{j\omega})=\frac{1}{2j}\left[e^{-j(N-1)(\omega-\omega_0)/2}\frac{\sin\left(N(\omega-\omega_0)/2\right)}{\sin((\omega-\omega_0)/2)}-\\\quad e^{-j(N-1)(\omega+\omega_0)/2}\frac{\sin\left(N(\omega+\omega_0)/2\right)}{\sin((\omega+\omega_0)/2)}\right]\tag{6}$$

I don't think that there is an analytical formula for the location of the maximum of the magnitude of $(6)$, but it is clear that due to the combination of the two terms in $(6)$, the maximum shifts away from $\omega_0$. That shift can become especially large for frequencies either close to DC or close to Nyquist. It is smallest for frequencies around half the Nyquist frequency (i.e., $f_s/4$). For small frequencies $\omega_0$, the actual frequency where the maximum occurs is smaller than $\omega_0$, whereas for large frequencies it is greater than $\omega_0$.


Information theory vs system theory ( commented in Determining linearity, causality, memory, and time invariability from a picture of the signal)

Yes that's most probably a typo; it's the system theory which deals with the concept of linearity, causality, etc.

Linear system theory deals with characterisations of systems in terms of its input-output relationship through applied signals. It describes how to compute the output in time and frequency domains as well as characterises systems according to their properties such as linearity, stability, causality etc.

Information theory deals with quantification of information (called as the entropy) contained within a source according to either its observed output (via statistical computation) or to its descriptive model (via probabiitistic computation).


There seem to be a number of hits in the literature from other domains that use various methods to segment the image (areas likely to belong to a filament and 'background') followed by filament identification  / extraction which you might want to have a look at anyway.

I was intrigued by the choice of PCA as a first indicator of extracting the local gradient and so this is a response with the "low hanging fruit" you can perhaps try (and in fact extremely easily) before jumping to the big guns of pattern recognition.

  ...the big question here is:
  What is the best way to find the filaments and/or their orientation at a given point?

Unless you have prior information (or assumptions) about the direction of the flow / filaments, it is impossible to derive their orientation. In other words, you can get some indication about the local orientation of (something that looks like) an edge, but not the orientation of the flow of the filament. You can assume, for example, that there is a particular intensity profile along the body of the filament which could then be used to determine "which part is the 'head' ".

The first thing to do (anyway) is to obtain the image gradient. The image gradient will give you the local slope of the intensity field (a.k.a image). All code here in Octave but only for demonstration. The operators are easily transferable across to scipy.

%Load the image
Q = imread('someImage.jpg');
%Resize it to 1/4 of its original scale, this is not essential but it helps with visualisation.
Qs = imresize(Q,0.25)
%
% Get the gradient
% [gradMag, gradDir] = imgradient(Qs);
% gradMag is the Magnitude, gradDir is the direction in degrees
% Break down gradDir to its components to process and visualise it
X = cos((gradDir/180.0).*pi);
Y = sin((gradDir/180.0).*pi);
%Let's have a look at the image with its gradient superimposed
imshow(Q2);hold on;quiver(X,Y);

The gradient of the (rescaled) image looks like this:

And if we zoom in (to the plot, not the image) around a bundle, we can see something like this:

So, predictably enough, the arrows are pointing away from the top of a ridge that signifies a filament and this can be used, to an extent, to derive its local orientation.

BUT, the bundle's intensity is not constant and the local gradient is not informed by any prior so that it "knows" which way to point.

Now, once you have moved from a scalar (Qs) to a vector field (X,Y), you might find other vector field operators useful as well, depending on what you are trying to do.

Can we do better?

Yes, we can apply a transformation to the intensity values which dampens the soft values and boosts the strong values. This will 'intensify' the ridges and basically result in a more well defined gradient later on. This is a modification of the image's dynamic range and is also known as 'expanding'.

An exponent is a good choice, it effects the following mapping:

Where we simply raise the pixel value to a power and take care of normalisation so that it does not overshoot the typical image limits.

This transforms the original (scaled) image to something like:

And conversely the local gradient now will be more well defined.

Now, if you are trying to detect broader structures which you would then infer (with some certainty) that they are formed by or belong to filaments, then you can do a couple of things. First of all, you could use Curl to identify "turbulent" (or absolutely no turbulence) areas and from there try to make sense of 'circularity' and then bundle together those areas.

The other thing you can do is something along the lines of line or circle detection through the Hough Transform. The Hough Transform is a "variant" of the Radon Transform. The latter is transforming an image from a $x,y$ plane to a $r, \theta$ plane by obtaining line integrals (the magnitude of which is $r$) at angle $\theta$ around the image. This looks like this (on Qs):

[R,xp] = radon(Qs);
imshow(R./max(max(R)));

The Radon transform converts lines (in Qs) to points. A line will produce a high value integral when $\theta$ is parallel to the bearing of the line. Conversely, arcs (and circles) are converted to arcs in the $r, \theta$ plane and if you follow intense arcs along the $\theta$ direction (corresponding to $x$ in a typical image) you will be able to identify the curled ridges that are formed by the filaments.

This is used by the Hough Transform which essentially operates on a thresholded version of the initial image (e.g. Qs) to detect lines and circles by looking for points with specific spatial relationships in the $r, \theta$ plane.

Hope this helps.


phase is the integral (w.r.t. time) of frequency.  and they are definitely not two independent processes.  back in, i dunno, the 30s, Armstrong implemented FM for the first time using AM and a synchronized carrier that was 90° out of phase with the AM.  so he was first accomplishing phase modulation (PM) out of AM (for which they had the existing technology - they didn't really have VCOs back in that day), then he implemented FM out of the PM by integrating the audio modulation signal and applying that to the PM modulation input.

this is all in a basic communications systems textbook like that of AB Carlson.


Not all LFSRs have the property of maximum cycle length. In fact, those who do have this property are based on a primitive polynomial of order n. It is a basic property of the extension field $GF(2^n)$ which is "spanned" by this primitive polynomial. The LFSR merely iterates through all of the powers of the polynomial primitive root.

Open a book on algebraic structures and cite the relevant part.


i have no idea what the "reconstruction fidelity term" is or what it's about.

Hermitian symmetry is a term usually applied to some form the Fourier Transform of a signal that is purely real.

for continuous-time, continuous-frequency Fourier Transform:

$$ X(f) = \int\limits_{-\infty}^{\infty} x(t) \, e^{-j 2 \pi f t} \ dt $$

if $x(t)$ is purely real (that is $\Im\{x(t)\}=0$ for all real $t$), then we know that there is this symmetry about $f=0$:

$$ X(-f) = X(f)^* = \operatorname{conj}\{X(f)\} \qquad \forall f \in \mathbb{R} $$

or

$$\begin{align}
\Re\{X(-f)\} &= \Re\{X(f)\} \\
\Im\{X(-f)\} &= -\Im\{X(f)\} \\
|X(-f)| &= |X(f)| \\
\arg\{X(-f)\} &= -\arg\{X(f)\} \\
\end{align}$$

Similarly, for the Discrete Fourier Transform:

$$ X[k] = \sum\limits_{n=0}^{N-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} $$

if $x[n]$ is purely real (that is $\Im\{x[n]\}=0$ for integer $n \in [0..N-1]$), then we know that there is this symmetry about $k=\tfrac{N}{2}$:

$$ X[N-k] = X[k]^* = \operatorname{conj}\{X[k]\} \qquad 1 \le k \le N-1 $$

or

$$\begin{align}
\Re\{X[N-k]\} &= \Re\{X[k]\} \\
\Im\{X[N-k]\} &= -\Im\{X[k]\} \\
|X[N-k]| &= |X[k]| \\
\arg\{X[N-k]\} &= -\arg\{X[k]\} \\
\end{align}$$

Now in my religion regarding the DFT, i insist that

$$\begin{align}
x[n+N] &= x[n] \qquad & \forall n \in \mathbb{Z} \\
X[k+N] &= X[k] \qquad & \forall k \in \mathbb{Z} \\
\end{align}$$

this is always the case (the DFT maps a discrete and periodic sequence of period $N$ in one domain to another discrete and periodic sequence of period $N$ in the reciprocal domain).  then the Hermitian symmetry takes a simpler form:

if $x[n]$ is purely real (that is $\Im\{x[n]\}=0$ for all integer $n$), then we know that there is this symmetry about $k=0$:

$$ X[-k] = X[k]^* = \operatorname{conj}\{X[k]\} \qquad \forall k \in \mathbb{Z} $$

or

$$\begin{align}
\Re\{X[-k]\} &= \Re\{X[k]\} \\
\Im\{X[-k]\} &= -\Im\{X[k]\} \\
|X[-k]| &= |X[k]| \\
\arg\{X[-k]\} &= -\arg\{X[k]\} \\
\end{align}$$

likewise you can identify Hermitian symmetry for the DTFT and relate it to the real-ness of the input to the DTFT. (do i have to do that one, too?)


Laplacian Operator with and without Diagonal Direction Elements in the Kernel

All of those are Discrete Approximation of the operator - Laplace Operator.
You can discretize it in any logical manner.
In the cases above, Istotropic means if you rotate it it looks the same.
The diagonal elements means the filter would be sensitive to changes which are in 45, 135, 225 and 315 degrees, it it means it will also amplify more noise.


How best to digitize a noisy phone recording for filtering?

Unless you have tight space constraints, I would stick to lossless formats like good old WAV. Phones do not pass much high-frequency content so you should be fine with a 8KHz Mono WAVE file. You will however want to keep the bit depth high, say to 24 bits, since you want to process the audio. Once you are finished with your processing you can compress with a lossy codec, such as OGG.

All codecs (coders/decoders) may be classified as lossy or lossless. Lossless codecs, as their name implies, allow for perfect reconstruction of the original signal. Lossy codecs, by contrast, merely approximate the source. By forgoing perfect fidelity, lossy codecs enable higher compression ratios, commensurate with the amount of degradation one is prepared to accept.

Perhaps related

A non-lossy aka lossless codec without patent-restrictions for audio?


How can I distinguish between two similar EEG signals?

yet another thing to try could be wavelet transform. Wavelets are powerful tools for analysis of signal. You could take a DWT of your signal and compare may be the most significant coefficients. Anyway but it would be helpful if you could post some sample data.


Understanding the Windowing Method in PSD Calculation

I don't really understand what do you mean by multiply them in the time domain and multiply them with window function. I think that you are trying to implement the Welch's PSD calculation. If so, steps should be:

Split your data into possibly overlapping segments of length $N$
Preferably remove the mean of each segment.
Window each of the segments by multiplying with window function $w$ of your choice.
Compute the Fourier Spectrum of such segment (you just calculated the modified periodogram)
Continue until all segments are processed, then average all of them. This will produce your PSD.
The last step is to scale your average by a proper factor.

For more details I am suggesting you refer to this document (Chapter 12 should be interesting for you, as well as previous ones):

  Heinzel G., et al. - Spectrum and spectral density estimation by
  the DFT, including a comprehensive list of window functions and some
  new flat-top windows

Possibly you can compare your results with pwelch function. For more theory, you can refer here (guys from Mathworks were most likely using this book). If you are using MATLAB, then this code should be helpful - very quick and dirty comparison I used years ago for pwelch. No overlapping, just two segments, but I am sure it could be a good starting point for you.

% Generate the random signal
N = 2048;
x = randn(2*N,1);
% Create the window and precompute scale
win   = hann(N);
s_win = norm(win,2)^2;
% Calculate the Welch's PSD
Pxx_welch = pwelch(x,win,0,N);
% Initialize the variable to store all segments
Pxx = zeros(N,2);
for k = 0:1
    x1 = x(1+k*N:(k+1)*N); % Take a segment
    xw = x1.*win;          % Multiply by window
    % Square the magnitude and scale
    Pxx(:,k+1) = abs(fft(xw)).^2./s_win;
end
% Average all segments
Pxx = mean(Pxx,2);
% Normalize accordingly to pwelch
Pxx = Pxx./(2*pi);
% Take only first half of frequencies
Pxx = Pxx(1:length(Pxx)/2+1);
% Preserve the energy (we just chopped-off half of spectrum), but do not
% modify first frequency bin as it is unique (mean value of a signal)
Pxx(2:end-1) = 2*Pxx(2:end-1);
% Plot the difference between two PSD's
plot(abs(Pxx_welch-Pxx))
title('Difference between pwelch and implementation')

Below is the plot of the difference between pwelch and the above implementation. One might notice that it is purely related to a round-off error.


Determining asymptotic stability using transfer function?

Poles in the right half-plane only cause instability if the system is a continuos-time system, because in that case, a pole $s_{\infty}$ contributes a signal

$$k\cdot e^{s_{\infty}t}u(t)\tag{1}$$

to the output signal, with some constant $k$, and with $u(t)$ denoting the unit step function. For $(1)$ to decay we require $\text{Re}\{s_{\infty}\}<0$, i.e., for a (causal) continuous-time system to be stable, all its poles must lie in the left half-plane.

In discrete time, however, a pole $z_{\infty}$ contributes a sequence

$$k\cdot z_{\infty}^nu[n]\tag{2}$$

to the output signal. Consequently, in order for $(2)$ to decay, we require $|z_{\infty}|<1$, i.e., for a (causal and) stable discrete-time system, all poles must be inside the unit circle.

Concerning the equivalence of poles and eigenvalues of the system matrix have a look at this answer.


Scikit-Image, Numpy, and Selecting Colors (python)

You can look at your picture in a different color space, e.g. the HSV space. In this space, each pixel has 3 components: hue, saturation and value. Hue defines the color's tone, saturation the color's strength and V the brightness of the pixel. For your particular problem, you can search for pixels that have e.g. a decent saturation (i.e. are reasonably colored and not just gray). The following code snippet uses OpenCV, but I believe you will be able to translate into scikit-image:

I = cv2.imread("/tmp/XApdz.jpg")  # Load the image

hsv = cv2.cvtColor(I, cv2.COLOR_BGR2HSV);  # Convert to HSV space
h, s, v = cv2.split(hsv);  # extract the single components

plt.figure(1, figsize=(10, 10))
plt.subplot(3,2,1); plt.imshow(h, cmap='gray'); plt.title('Hue')
plt.subplot(3,2,2); plt.imshow(s, cmap='gray'); plt.title('Saturation')
plt.subplot(3,2,3); plt.imshow(v, cmap='gray'); plt.title('Value')

mask = (s > 80).astype(np.uint8);  # Thresholding in the Saturation-channel
plt.subplot(3,2,4); plt.imshow(mask); plt.title('mask')

disk = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(3,3))  # remove very small regions
opened = cv2.morphologyEx(mask, cv2.MORPH_OPEN, disk);
plt.subplot(3,2,5); plt.imshow(opened); plt.title('Opened mask')

I2 = I.copy()
I2[opened.astype(bool), :] = 0; # Set the pixels to zero, where
plt.subplot(3,2,6); plt.imshow(I2)


How do you calculate/compute/determine the spectral coverage of a frequency spectrum?

Spectral coverage, I assume, means the Bandwidth required to transmit the incoming signal. A simple formula(approximate) for the Spectral Efficiency of Memoryless Modulation Schemes(MSK is not one) which gives the ratio between Data Rate and Bandwidth required is:

R/W = (2*logM)/N

where R is the incoming Data Rate and W is the required Bandwidth; log is in Base 2, M is the number of points on the Signal Constellation and N is the number of dimensions required to use the particular signalling.

FSK requires as many as Dimensions as the number of Points you want to transmit.
ASK is one-dimensional.
QPSK is two-dimensional.
MSK is more spectrally efficient than the other 3 and I think the formula for it's Power Spectral Density and its 3dB value can be found in all Digital Communication Texts.


Frequency conversion in IIR lowpass-to-lowpass transformation

The prototype filter is not uniquely defined by the two band edges of the filter to be designed. So you have to choose one of the two edge frequencies of the prototype filter. E.g., if you choose $\omega'_p$, you can compute the value of $\alpha$ according to
$$\alpha=\frac{\sin[(\omega'_p-\omega_p)/2]}{\sin[(\omega'_p+\omega_p)/2]}\tag{1}$$
where $\omega_p$ is the desired pass band edge. Now that you have $\alpha$, the stop band edge of the prototype filter can be computed from the desired stop band edge $\omega_s$:
$$\omega'_s=-\angle\left(\frac{e^{-j\omega_s}-\alpha}{1-\alpha e^{-j\omega_s}}\right)\tag{2}$$

EDIT: the lowpass-to-lowpass transformation
$$z^{-1}\longrightarrow\frac{z^{-1}-\alpha}{1-\alpha z^{-1}}\tag{3}$$
transforms a prototype edge frequency $\omega'_s$ in the following way:
$$e^{-j\omega'_s}=\frac{e^{-j\omega_s}-\alpha}{1-\alpha e^{-j\omega_s}}\tag{4}$$
From $(4)$ it follows that
$$-\omega'_s=\angle\left(\frac{e^{-j\omega_s}-\alpha}{1-\alpha e^{-j\omega_s}}\right)\tag{5}$$
which is equivalent to $(2)$.
Note that Eq. $(5)$ is different from Eq. $(8.71)$ in the book because $(8.71)$ is derived from a lowpass-to-highpass transformation.


Converting complex baseband signal to passband signal

Your signal $z$ is centered at 20 MHz, but because your sampling rate is 30.72 MHz, 20 MHz is beyond the Nyquist frequency of $f_s/2$ = 15.36 MHz and therefore aliases to 20-30.72 = -11.72 MHz.  The signal straddles $\pm f_s/2$, which is not what you want.  To center a 20 MHz wide signal at 20 MHz, you must sample at least $2*(20+20/2)=60$ MHz, or roughly twice your current sample rate.  Thus you want to interpolate the signal to a higher sample rate first.

The other issue is that your signal is still complex.  Maybe that is what you want, but typically if you convert to passband it is because you need a real signal - e.g., for D/A conversion.  In that case you also need to take the real component of the signal after the complex upconversion.

If you do need a passband signal for D/A conversion, then you might want to interpolate by more than 2 to make your D/A reconstruction filter easier - but that is another topic entirely.


Remove things you don't want

Since the camera is static, you might want to use a background remover first. I found that the standard one provided with OpenCV works pretty well. I create it like this in the Android OpenCV SDK (you can play with the parameters) :

backgroundSubtractor = new BackgroundSubtractorMOG(3, 4, 0.8);

Then, apply it to each image in your video (no need to re-instantiate it unless you move the camera, it will learn from the sequence of images).

Computing the motion

First, using the Hough Transform for fast-moving shapes isn't very useful. The Hough transform requires you to be able to detect the shape's boundary with some accuracy, which obviously you can't, thanks to the motion blur.

I would suggest looking into algorithms for computing the optical flow. I don't know how close EmguCV is to the original OpenCV, but you could try calcOpticalFlowPyrLK (computes the flow for a set of points in the image, usually obtained using the goodFeaturesToTrack method) or calcOpticalFlowFarneback (computes a dense flow for the whole image, but more CPU-hungry).

Additionally, there have been some people working on CUDA/GLSL versions of Lucas-Kanade (the first one) if you really want to have a fast, dense flow algorithm, but I haven't seen a lot of useful implementations.

See the OpenCV documentation here : http://opencv.willowgarage.com/documentation/cpp/motion_analysis_and_object_tracking.html#calcOpticalFlowPyrLK


Can you increase frequency resolution of FFT without increasing window size?

If you really insist on using FFT (rather than parametric methods, which wouldn't suffer from time/frequency trade-offs), you can fake a much better resolution by using the phase information to recover the instantaneous frequency for each FFT bin. Partials can then be detected by looking for plateaus in the function giving instantaneous frequency as a function of FFT bin index. The common implementation of this technique as described in this paper will "cost" you one extra STFT (instantaneous frequency is recovered by operations on the STFT of the signal, and STFT of the derivative of the signal).

See for example the ifgram function in this Matlab implementation of sinusoidal modeling of audio signals.

Note that this won't help resolving two partials falling into adjacent FFT bins. It will just provide a much more accurate frequency estimate than just converting into a frequency the FFT bin index of a spectral peak.


How to combine a perfect signal with a limited dynamic range with a poor one with high dynamic range?

In response to your question (1):

Thinking along the lines of your initial idea, there is a way to avoid regression and do the following approximate procedure instead:

Without loss of generality, assume that $\epsilon$ has zero mean (if not, then just absorb it into the constant $\beta_0$). Also assume it is i.i.d and has finite variance. Suppose you obtain two sets of data for $(f,g)$ pairs when $h=1$. I will denote these as $\{(f_{1,i}, g_{1,i})\}_{i=1}^N$ and $\{(f_{2,i}, g_{2,i})\}_{i=1}^N$.

Let $\bar{f_1} = \frac{1}{N} \sum_{i=1}^N f_{1,i}$ and similarly define $\bar{g}_1$, $\bar{f}_2$ and $\bar{g}_2$. If $N$ is sufficiently large, we can now claim that, at least approximately, the noise $\epsilon$ can be averaged out and we can write:
$$
\bar{g}_1 \approx \beta_0 + \beta_1 \bar{f_1}
$$
and
$$
\bar{g}_2 \approx \beta_0 + \beta_1 \bar{f_2}.
$$
Now we only need to solve two equations in two unknowns to get the scale parameter $\beta_1$ and offset $\beta_0$, provided $|\bar{f}_1 - \bar{f}_2|>0$.

Of course, whether this will work in your application will depend on whether the variance of $\epsilon$ is small enough (remember: averaging will reduce the variance by a factor of $1/N$), and whether the $f$'s are really noise free when $h=1$.

In response to question (3):

The same procedure should work in case $f$ is a quantized version of the true velocity. The argument here is that quantization introduces a small amount of noise that is approximately uniformly distributed with zero mean. It gets averaged out for large $N$, or so I hope.


How can the FFT be used for estimating linear chirp parameters?

Why use an FFT for this? I think there are better simpler ways to do this. For example, you can simply measure the distance between zero crossings to get a first order estimate for frequency vs time and then refine with reconstructing the sweep and looking at the difference of the original and the reconstructed signal.


A continuous-time signal $x(t)$ is sampled to create an ideal discrete time signal $x[n]$ by just taking the instantaneous value of the continuous-time signal at the sampling instants. If you sample the signal at an interval $T_s$, then:

$$
x[n] = x(nT_s) = [0, x(T_s), x(2T_s),\ \ldots \ ]
$$

In your case, $T_s = \frac{1}{8000}\text{sec}$. In practice, there are no ideal samplers, so this model doesn't strictly hold, but it is close enough for most applications and is the model that you would typically see used in introductory signal processing textbooks and courses.


How do I numerically calculate a function from its noisy gradient?

You can filter either the gradient itself or the result, $s$.  You would need to know the characteristics of the true gradients fairly well to know what the frequency bandwidth is.  At that point you could design a low-pass filter that would preserve the signal but get rid of higher frequency noise.


How can initial conditions be taken into account to calculate a system's terminal value using the final value theorem or some other technique

Ah ha!  There's a lot of obfuscation in the problem as stated, plus a bit of over-emphasis on the final value theorem.

$$H(z) = \frac{z^{-1} \phi / (1 - \phi z^{-1})}{1 - z^{-1} \psi / (1 - \psi z^{-1})} \tag 0$$

OK, so far so good.

It arises from the following feedback system
$$Y(z)(1-z^{-1}) = z^{-1}F_{\psi}\cdot Y(z) + z^{-1}F_{\phi}\cdot X(z) \tag{1}$$
where $F_{\alpha}(z)$ is a 1st order high pass filter with parameter $\alpha$ of the form
$$F_{\alpha}(z) = \frac{\alpha (1 - z^{-1})}{1 - \alpha z^{-1}} \tag{2}$$

This was bothering me yesterday, but I didn't notice the essential problem.  The general form of the filter described in (2) is first-order, and the delay operators ($z^{-1}$) on the right side of (1) potentially add a second pole.
The original transfer function (0) is in a screwy form, but it appears to be second-order.  So the whole problem manages to worm its way through the pole-zero cancellation check by virtue of being stated in such a non-standard way.  This question deserves to be an extra-credit problem in a control systems class, or problem 10 out of 10 in a final by a particularly sadistic teacher, or perhaps a routine problem if the instructor believes in homework sets that no one will be able to fully solve.
The problem is that there's a pole-zero cancellation that snuck right by me and everyone else.  It's evident in the left-hand side of (1), where the derivative of $y_k$ is the subject of the equation.
So the reason that you can't solve this problem as stated using the final value theorem is that you cannot adequately represent the system using a transfer function.  There may be some way to save this within transfer function notation, but I just tried and failed at the first step, so I'm going to do it in state-space.
So, first, shine some light on this obfuscated* problem.  I'll do that by substituting (2) into (1):
$$Y(z)(1-z^{-1}) =
     z^{-1}\frac{\psi(1-z^{-1})}{1 - \psi z^{-1}}Y(z) +
     z^{-1}\frac{\phi(1-z^{-1})}{1 - \phi z^{-1}}X(z) \tag 3$$
Now the pole-zero cancellation becomes painfully evident: we're expressing a filter of the form $1 - z^{-1}$ on the left-hand side, but when we divide through by $1 - z^{-1}$, it no longer appears in the equation.
Oh joy.
So, drawing this out as a feedback loop (and changing from $z^{-1}$ to $z$ notation, because I'm a control systems guy):

Notation is going to get a bit awkward here.  Bear with me.  The transfer function in (2) can be realized as a state-space system (where $x$ is a state variable, not the system input):
$$\begin{aligned}
x_k = \alpha x_{k-1} + (\alpha - 1)u_k \\
y_k = -x_k + u_k
\end{aligned}
\tag 4
$$
So (awkwardly) noting that $x_{1,k}$ is the first element of the state vector $\mathbf{x}$ at time $k$, but is not the system input $x_k$, then the 'top' input to the summing junction is (4) with $\alpha = \phi$ and an extra delay:
$$\begin{aligned}
\begin{bmatrix}x_{1,k} \\ x_{2,k}\end{bmatrix} =
\begin{bmatrix} \phi & 0 \\ 0 & 0\end{bmatrix}
  \begin{bmatrix}x_{1,k-1} \\ x_{2,k-1}\end{bmatrix} +
\begin{bmatrix}\phi - 1 \\ 1\end{bmatrix} x_k \\
u_{1,k} = x_{2, k-1} - x_{1,k-1}
\end{aligned}
\tag 5
$$
The output of the system, $y_k$, is itself just a state variable (oh joy -- not sarcastically this time).  So just let $$y_k = x_{3, k} \tag 6$$.
The filter in the feedback path is again described by (4), with $\alpha = \psi$ and an extra delay:
$$\begin{aligned}
x_{4,k} = \psi x_{4,k-1} + (\psi - 1)x_{3, k} \\
u_{2,k} = x_{4,k-1} - x_{3, k-1}
\end{aligned}
\tag 7
$$
Finally, the integrator in the forward portion of the loop is:
$$x_{2,k} = x_{2,k-1} + u_{1,k} - u_{2,k} \tag 8$$
This goes together into a fourth-order system:
$$\begin{aligned}
\mathbf{x}_k = \begin{bmatrix}
\phi & 0 & 0 & 0\\
0 & 0 & 0 & 0 \\
-1 & 1 & 0 & -1 \\
0 & 0 & \psi - 1 & \psi
\end{bmatrix}
\mathbf{x}_{k-1} + \begin{bmatrix}\phi - 1 \\ 1 \\ 0 \\ 0\end{bmatrix}x_k \\
y_k = \begin{bmatrix}0 & 0 & 1 & 0\end{bmatrix}
\end{aligned}
\tag 9
$$
Recast that into $\mathbf{x}_k = A \mathbf{x}_{k-1} + B\,u_k,\,y_k = C\mathbf{x}_k$, and you can test the $A$ and $C$ matrices to see that it is observable, at least for some combinations of $\phi$ and $\psi$ (it is, I checked).  Since there's pole-zero cancellation, we can assume that it isn't controllable (you can check that for yourself).
Because it's observable, we can consider the system
$$\mathbf{x}_k = A \mathbf{x}_{k-1} + \mathbf{x}_0 \delta(k)$$
The response of this system to any given $\mathbf{x}_0$ is $$Y(z) = C(\mathbf{I}z - A)^{-1}\mathbf{x}_0$$  Now you can find the system response for each initial value in $\mathbf{x}$ using the final value theorem.
* I looked it up -- "obfuscate" literally means "to darken, to put into shadow".  So shining a light on the thing really is the antonym of "obfuscate".


Ultrasound pulses and reflection problems

It would appear that your problem lends itself very nicely to using a CDMA scheme.

Let starts with some properties of (DSSS) CDMA. (Direct Sequence Spread Spectrum, Code-Division Multiple Access). Its a mouthful, but it is really easy to implement.

In CDMA, your pulse (at baseband) is actually made up of many concatenated 'chips' as they are called. The chips are just 1s or -1s, of a fixed duration. For example, your chipping sequence might be [1 -1 1 -1 -1 -1 1]. You would use this chipping sequence to modulate your carrier.

However, you cannot just make up your chipping code. What you want to do is use chipping codes that have the very nice property, that their autocorrelation function is a delta function like so:

(Equivalently, their power spectral density is white). For example, you can look into using Barker Sequences as your chipping code, (usually used in radar), or you can also look at using Gold Codes. Practically speaking however, this means that you get the maximum correlation score in your receiver, ONLY when the receivers' code, exactly lines up with the transmitted code, and zero otherwise.

How does this help you? In your receiver, you would be running a correlator continuously. The correlator would be performing a running dot-product of its own local code, with whatever is received. Now imagine that you receive a transmitted waveform from your pen, and a second waveform from a reflection. As your receivers' correlator runs, it will give a peak when its own codeword exactly alligns with your code from the pen. This will cause your detector to 'lock' onto that specific delay value. Now, here is where you reap the benefits of a near-delta autocorrelation function of your code: The reflected signal will also be present, and will also have its dot product taken with the receivers' locked code, but it will give zero, or near zero score, since it is orthoginal or near-orthogonal to the delayed code that your receiver has already locked onto.

In contrast, if you had send out a un-coded carrier pulse, you would be at the mercy of constructive or destructive interference throwing off when exactly your pulse peaked at the detector level of your receiver, and thus get erroneous TDOAs.


Sampling Rate of the daq tool tells you how many measurement samples will be collected per second. So, if the sampling rate $f_s=1Msps$, it means that your daq tool is collecting $10^6$ samples per second.


Signal equation using signal waveform and Fourier series?

Did you check the usual equations?

Responding to:

  The steps I should follow are

  Find $f(x)$ from the signal waveform
  Compute $a_0$, $a_n$, $b_n$
  Put above values into the (2) and see whether I get (1) or not.

Essentially, yes.

For your equation to relate to the $f(x)$ in the image above you may want to think about
$$
w_m = \frac{2\pi 5 m}{100} = \frac{\pi m}{10}
$$

Also note that there needs to be an integer (or rational) ratio between $f_1$ and $f_2$ in your equation (1).


Transforming digitized noisy signal before applying cross-correlation

To offset or not is simply in your definitions on how you want to do the math involved and can be convenient for further processing.  Subtracting a constant, or otherwise scaling the waveform, does not change the signal to noise ratio.

Correlation is to multiply and accumulate, and the cross-correlation and auto-correlation functions show this correlation as the sequences are shifted in time relative to each other.  Cross-correlation is between two different functions, and auto-correlation is between a function and itself.

The cross-correlation of x[n] and y[n], with a shift index m, can be given as:

$$\rho_{xy}[m] = \Sigma_n x[n]y[n+m]$$

(Note: For complex signals, the product must be complex conjugate product)

With that said, consider this very simple example using an 11 chip barker code to demonstrate my first statement. This Barker code has a similar property to C/A codes in that a correlation with shifted versions with itself (autocorrelation) will have a very low correlation relative to when the code is aligned.

1 0 1 1 0 1 1 1 0 0 0

This code written with 1's and 0's has a mean of approximately 1/2, so you can optionally remove the mean and similarly for convenience scale by a factor of two by mapping 0 to -1 such that the code is:

1 -1 1 1 -1 1 1 1 -1 -1 1

This results in the following correlation when the sequences are aligned in time (which could be used in direct-sequence-spread-spectrum, DSSS, to send a data symbol with value = "1"):

And the following if the code was inverted (which could be used in DSSS to send a data symbol with value = "0"):

And the following for any other shift (I did a rotational shift of one sample here, corresponding to m = -1):

With the following autocorrelation result as a plot:

We could similarly instead perform the correlation with the following operations for multiplication (basically XOR function, could do the same with XNOR and get an inverterted result):

1 x 0 = 0

0 x 1 = 1

1 x 1 = 1

0 x 0 = 0

Resulting in a similar result with an offset and scaling in the correlation output instead of an offset and scaling the correlation input. So in this case the maximum correlation would be 11, the maximum correlation when inverted is 0 and the autocorrelation for any other offset would be 5. Double this and subtract 11 and you get the same result as before.

(Note unlike maximum-length pseudo-random sequences which will also cross-correlate to -1 at all offsets when you use +1 and -1 values, C/A codes have cross-correlation values of either 1, +63 or -65, which is still significantly less than the 1023 value when aligned).

The reason this works for GPS in the case of low SNR is when you add each of the 1023 chips with noise present, if the noise on each chip is independent ("white noise"), then the signal magnitude will accumulate by a factor of $1023$, but the standard deviation of the noise will only accumulate by a factor of $\sqrt{1023}$ thus in terms of SNR you get a $10Log_{10}(N)$ processing gain. Consider a simpler case of adding two independent Gaussian random variables with equal distribution and equal but non-zero mean: The mean will double but the standard deviation will only increase by the $\sqrt{2}$.


what does "derivative" means in image processing?

Look at the numbers in the filter kernel in just 1 dimension for simplicity. For a Sobel and  Prewitt matrix you have something that roughly looks like this

[-1,0,1].

Convolving this with your image basically computes the difference between the pixel values of the neighboring pixels. You apply 0 to the current pixel, 1 to the pixel on the right and -1 to the pixel on the left. This gives a first order difference:

next pixel - previous pixel,

i.e. first derivative.

But now look at a Laplacian operator. It looks something like [1, -2, 1]. This computes the difference of differences. To see how, note that

[1,-2,1] corresponds to

next - 2 x current + previous

i.e.

next - current - current + previous

i.e.

(next-current) - (current-previous)

Now notice how this is a diference of differences. (next - current) is like a 1st derivative. (current - previous) is like 1st derivative. Their difference is like a 2nd derivative.


All these notation should be considered within the context of the development.

So $r_d[k]$ may refer to the auto-correlations sequence of the desired signal whereas the $r_{dx}[k]$ is the cross correlation between the input $x[n]$ and desired signal $d[n]$. So they are two different things.

The latter is used in Wiener optimal filtering to solve the Wiener-Hopf equations (aka The Normal Equations) but I cannot recall where the former is used, can you cite the formula, chapter etc. that makes use of it?


How to draw a vector between light and dark spots on an image

The keyword here is blob detection. Search OpenCV blob detection and I'm sure there are many examples.

One of the common methods is to use the Laplacian of Gaussian operator.
The Laplacian is a measure of the strength of the second-order derivatives. If you think of the blobs as peaks and valleys on a 3D surface, at the middle of the blobs the slope (derivative) of the surface is changing (second derivative) rapidly from negative to positive and vice-versa.

The Gaussian is there to average out this measurement over a wider area. Any area averaging kernel could be used, but the Gaussian is special* in that increasing the size of the kernel does not introduce any new features. Using multiple kernels creates a scale space, and allows you do scale invariant blob detection.

*I think the Poisson kernel can create a scale-space too with a slight rule relaxation.


The topic is to the best of my knowledge an area of open research. The research contribution which helped me most in solving my problem was: The spectral analysis of cyclo-non-stationary signals. Specifically, the key seems to be:

"In this paper, it is shown that the optimal way of describing
cyclo-non-stationary signals is jointly in the time and the angular
domains."


Quasi-static channels can be said as "block-wise" time (in)variant. For example, your channel has no variations over time (not in delay domain) for some time-period, e.g., 1 ms, but the channel may change after 1 ms.


No; the system given by

$$ y[n] = x[n] + n $$ is time-varying, due to the added term $n$.

Your mistake is in the line :
$$\ y_2[n] = x_2[n] + n = x[n-k]+(n-k) $$ which should be instead

$$\ y_2[n] = x_2[n] + n = x[n-k]+n$$ and therefore implies that

$$y_2[n] \neq y[n-k] $$.


Deriving time-scaling property for Fourier Series

It looks more like zero padding to create up-sampling resolution to me.

Prove
x[n/m+N]=x[n/m], since your argument is now n/m and not n.


Low-pass filtering bounces when signal doesn't start at 0

The problem is due to the fact that initially, the filter memory is filled with zeros (if not otherwise stated). Since the filter output is a weighted average of the filter memory and the input, it will start its output near zero, when the initial state is set to zero.

As Atul suggested, you can overcome this by padding 30 to the beginning of the sequence (to fill the memory before you actually filter the sequence of interest). However, the amount of padding depends on the filter and how fast it adjusts to the steady state (depending on filter order and e.g. Q-factor)

As a more elegant solution, you can initialize the filter state such that its first sample already outputs 30. You can calculate the filter state as the steady state response as to a unit step input. E.g. python has this functionality already built-in:

signal = np.linspace(30, -30, 1000)
sigma = 1
signal = signal + sigma * np.random.randn(len(signal))

plt.plot(signal)
b, a = scipy.signal.butter(5, 0.1)
zi = scipy.signal.lfilter_zi(b, a)

filtered = scipy.signal.lfilter(b, a, signal, zi=30*zi)[0]
plt.plot(filtered, lw=2)


There are a number of things that could be going wrong; you haven't provided enough information to definitively diagnose your the problem. (A link to the wav file you obtained would be useful, so we can try it out ourselves.) Here are some observations:

You're assuming the file is uncompressed 16-bit little-endian samples, and ignoring the presence of the header. This could be wrong, either because the format is not what you think, or the header is not a multiple of 4 bytes long. Instead, use an audio-editing program such as Audacity to read the file and write it out again as a raw file where you have specified the parameters. That, or use a WAV-reading library in your program.
Your lengths don't match up. If you read 1024 bytes, then you have 512 shorts, or 256 complex-shorts — but you say your FFT size is 512, not 256. Either your program is not as you say or it might be reading uninitialized memory.
Your output appears to have four almost-copies of the spectrum. Two mirrored copies are explainable if you are misreading the input such that the I and Q components are not aligned properly — a symmetric spectrum is to be expected if one of I or Q is zero, constant, or otherwise not properly aligned with the input.

The other copy might be due to an error in your FFT shift code or graphics code, or because you are trying to display both components of the complex FFT output separately rather than taking the magnitude.
I would suggest adjusting the scaling of your color shading so that it is more similar to the HDSDR output, and increasing the FFT size (bin count). (You may also need to take the log of the magnitude.) This will make the result look more like HDSDR's, which will help you match up parts of the spectrum.

  Do I need to process the IQ signal before somehow (demodulate it ?)

No. The sorts of processing that could be done to refine a panadapter display would be done after the FFT. The only things to do before would be things like IQ imbalance correction, which are just as appropriate for a signal being demodulated and played as audio (or whatever sort of information) as one going to the panadapter.


Well my friend, if you read the following B&K tutorial about the Order Tracking, then you should encounter so called Vold-Kalman Order Filtering. Personally I don't like it (rather prefer Chirp-Fourier based methods for extracting the response of the system).

You can find plenty implementations in the internet. When you look for MATLAB, then please refer to:

Vold-Kalman order tracking code

Second generation Vold-Kalman Order Filtering.

And most certainly I suggest you to search the bottom of the prof. Jiri Tuma website: click!

For more theory please take a look at - I'm sure you will like it:

Characteristics of the Vold-Kalman Order Tracking Filter

Vold-Kalman Order Tracking Filtration

Main Principles and Limitations of the Current Order Tracking Methods

Realization of the Vold-Kalman Tracking Filter - A Least Squares problem

Advanced Vold-Kalman Filtering Order Tracking - do not be scared, it's not in chineese ;)


In fact both plots are the same.

Your freqz plot uses a linear discrete frequency, $\omega$ axis scaled between $0$ and $1$ with units of $\pi$ rad per sample. For 200 Hz sampling frequency, a 60 Hz notch would be shown at $w_c = 0.6 \pi$ rads, which is put at 0.6.

The bode plot, however, uses a logarithmic scale for the continuous-time radian frequency $\Omega$. Given that you have a $f_c=60$ Hz of notch frequency, its radian value is $\Omega_c = 2\pi \times  60 \approx 377$ radians per second.

So both plots are showing the notch at the correct place.


Discrete-time signals are generally written in the time domain as $x[n]$ where $n$ is an integer.

We want to be able to see how discrete-time systems can be used to modify these signals:
$$
y[n] = \sum_{m=1}^M a_k y[n-m] + \sum_{p=0}^P b_p x[n-p]
$$

The above equation is a linear, constant coefficient difference equation.

Such equations can be hard to manipulate for a specific $x$ to figure out what $y$ is.

The $z$ operator is used to attempt to simplify the analysis of such systems. $z$ is the forward-shift operator: it shifts a signal forward by one time instant. $z^{-1}$ is the backward-shift operator: it shifts a signal backward in time by one time instant.

Using this operator, we can rewrite the equation as:
$$
Y(z) = Y(z) \sum_{m=1}^M a_k z^{-m}  + X(z) \sum_{p=0}^P b_p z^{-p}
$$
which, upon rearrangement, we can get:

$$
Y(z)/X(z) = \frac{\displaystyle\sum_{p=0}^P b_p z^{-p}}{\displaystyle1 -  \sum_{m=1}^M a_k z^{-m}}
$$

which is now a discrete-time system that we can analyze in the $z$ domain, rather than the time domain.

One thing this buys us is the ability to figure out what happens when many discrete-time systems are put together: we just multiply the discrete-time system descriptions (transfer functions).  Doing the same thing in the time domain is doable, but very tedious.


What is averaging and how can it be done?

In general, EEG, fMRI (and also MEG, SPECT and PET datasets -the so called "functional modalities") are obtained from a subject (e.g. a human being) while it is engaged in one or more "activities".

This "activity" could even be "Stay at rest with your eyes closed" or "try to solve this equation" or "simply look at the images on this screen".

Usually, more than one subjects will be engaged in a study and these subjects may be divided into groups according to some criterion to compare and contrast their responses (e.g. "Controls" (subjects that are "normal" or do not suffer from a disease) and "Diseased" (subjects that have been verified to suffer from some disease)).

Both of these modalities, in addition to sampling time they also sample space AND each fMRI / EEG dataset is itself a sample from a population of subjects mentioned above.

Having described this, you can almost see the dimensions across which you would apply averaging and why:

You can do (the common) temporal averaging (like a moving average or median filter) over each channel (for the EEG) or voxel time series (for fMRI) to reduce the effect of noise.
You can do Spatial averaging in a number of different ways and for different purposes but all that spatial averaging entails is summing together the responses from $N$ individual channels ($N$ individual vectors) and diving them by $N$ to get a mean response.

For example, perhaps you would like to obtain the average response corresponding to some brain region (occipital, parietal, frontal, other). Then you can select the $N$ vectors (channels) corresponding to those regions and obtain ONE vector by averaging them (for example, if your channels are in X(channelN, timeSampleM) you can obtain ONE average by doing mean(X) which will return one vector of size 1 x timeSampleM with the mean of each column (that corresponds to the sample taken by each channel at some time instant $M$).

Another thing you might want to do, is obtain the average response of ONE subject depending on the activity they were carrying out. For example, you would like to see the differences in the average response of a subject's brain between the conditions of REST and LISTENINGTOMUSIC. In this case, you need to know which datasets correspond to which activity and you simply obtain the average of those $X_{k,l}$ (where $k$ denotes the dataset and $l$ denotes the activity).
You can do spatial averaging over a population (and / or "activity" they were carrying out). So for example, you might want to see what is the average response of the brain at condition REST between two subject groups (Group A, Group B). In this case you select all datasets from all subjects in Group A and obtain an average dataset and similarly for Group B. You can do the same (but changing the membership) to obtain the average response of ONE group between different conditions.

(And of course, you could apply any combination of spatial / temporal / group / condition averaging depending on what you are trying to observe)

Finally, another area where averaging is inherent to the data processing is that of Evoked Potentials.

For example, imagine that you would like to examine the response of the brain to a burst of audible noise. In this case, you provide a stimulus to the subject and record its response for a number of times. After that, you average all recordings (alligned at the stimulus instance) and produce ONE signal that represents the response of the brain. A famous such signal is the P300.

As you can see, although the operation of averaging is the same (from a mathematical point of view) there are many different ways to apply it depending on the objective. Therefore, I would recommend going over a specific paper carefully to understand

a) why they apply averaging
and b) how is averaging applied.

Furthermore, i would like to recommend to you the following resources that you might find useful and which cotain far more detailed information about EEG and fMRI analysis (incudling averaging) in general:

Electroencephalograhy by Lopes Da Silva (obviously for EEG, key
chapters to focus on are 2,4,14,53,58,59)
Statistical Parametric Mapping by Friston et al (for a very
thorough treatment of the analysis of fMRI images. Earlier versions
of that book (and of the also excellent "Human Brain Function" are
available for free from here (you need to have a close look at
chapters 3-5 and especially 5 to understand the experimental side))
Scholarpedia articles on Neuroimaging (This is a very good
resource in general)


What is "Equivalent Time Sampling" and what is it good/used for?

Note that the Nyquist criterion applies to bandlimited signals; but the so-limited band does not need to be baseband, or even contiguous, just not overlapping (aliased) with any other sub-bands after folding.

Periodic sampling fan-folds an infinite spectrum in to "fans" with a width of half the sampling rate.  Any frequency bands, after this fan-folding, that do not overlap with any other frequency bands in the spectrum when so folded can be reconstructed, according to one form of Nyquist theorem.

In practice, any random phase jitter in undersampling will cause the bands too far above Fs/2 to become more phase noise than signal.  If the phase noise is mostly white and uncorrelated with the signal, then sampling for a much longer time frame may allow averaging out this phase noise sufficient to reconstruct a repetitive narrow-band signal, with some probability of accuracy for a given extended sampling time window.

So, some oscilloscopes have a mode where they undersample a bandpass-filtered narrowband input, then store and composite a much longer sampling time by folding the samples down to baseband.

In addition, signal reconstruction of a bandlimited signal may be possible with non-periodic sampling of known sample times, as long as the local density of samples is high enough.  This is numerically much less stable, as estimating a reconstruction of the spectrum or waveform from a minimum number randomly spaced samples becomes more of a sensitive linear optimization problem, rather than a closed form transform or convolution.  However some oscilloscopes might use a pseudo-random non-periodic undersampling, but averaged over a long time period, and let the user determine when the density of samples "looks good enough" when reconstructing an image of an undersampled repetitive waveform.

Randomized sampling can avoid certain aliasing problems, such as perfectly periodic sampling not allowing reconstruction of an arbitrary phase sinusoid with a frequency exactly at Fs/2 (it's smooshed at the creases in the fan-fold).


Binary and Quadrature Phase Shift Keying

(1) The phases are indeed relative. The numbers 45, 135, 225 and 315 are relative to the transmitter's oscillator. The receiver, of course, has a different oscillator, so it will need to do some procesing to estimate the phases.

Related to this subject, there is also differential modulation, where 1s and 0s are transmitting by changing the phase relative to the previous symbol.

(2) This identity is useful to understand I/Q modulation:
$$A_I \cos(2\pi f t) + A_Q \sin(2\pi f_c t) = S_{IQ} \cos(2\pi f_c t - \theta_{IQ}) $$
where
$$ S_{IQ} = \sqrt{A_I^2 + A_Q^2}$$ and $$\theta_{IQ} = \tan^{-1}\frac{A_Q}{A_I}.$$
So, as you can see, you can achieve any phase $\theta_{IQ}$ by adding the I and Q signals with appropriate amplitudes.

(3) Binary PSK is just an I signal. Since the two phases in BPSK are zero and 180, the Q portion of the signal is always zero.


For time-based systems, I understand that it is difficult to imagine a memory of the future. But for general systems, $-t$ and $t$ are just left and right (think of a spatial system). Other discussions are in LTI system $y(t)=x(t−T)$ with or without memory, What is a memory less system?, or A question about the concept of the time.
By definition of invertibility, $\mathcal{L}^{-1}$ is such that  $\mathcal{L}^{-1}( \mathcal{L}(x))=x$. But also that $\mathcal{L}( \mathcal{L}^{-1}(x))=x$ (by the way, derivatives and integrals are not inverses). Let us suppose the converse: $ \mathcal{L}^{-1}$ has no memory. Hence $\mathcal{L}^{-1}(x[n]))$ can only use the present state, and $\mathcal{L}$ as well to yield $(x[n])$.
So, if $\mathcal{L}^{-1}$ is memoryless, $\mathcal{L}$ is memoryless as well. By contraposition, the converve is true


MATLAB gradient derivative troubleshooting

You can calculate the derivative using the diff function and then dividing by the time interval between the "samples" in your vector like this:

dA = diff(A) / 286

The diff() function calculates the difference between an element and the previous one in the vector. With this you'll be calculating the derivative for each element in the vector as the increment in A divided by the increment in x. Check out the help for the diff() for more info.


Why a pair of complex-conjugate zeros provide a nulling filter? (FIR filter case)

The unit circle on the z-plane represents the frequency axis, similar to the imaginary axis $j\Omega$ on the s-plane for the Laplace Transform in the continuous time case. So the frequency response of the system is given by $H(z)$ when $z= e^{j\omega}$ with $\omega$ going from $0$ to $2\pi$ representing the normalized fractional radian frequency (which is the continuous time radian frequency $2\pi f$ divided by the sampling rate $f_s$.

That said, any zero on the unit circle will create a null in the frequency response. With the OP's case of complex conjugate zeros (resulting in a real response), two nulls would result as shown. The location of the zero if on the unit circle, is the fractional radian frequency where $H(z) = 0$, thus called a "zero".

If the zero is not on the unit circle, the null will not be zero but will be lower the closer the zero is to the unit circle for that frequency.

This may be clearer from the plot below showing the frequency response for a 2 point moving average filter, which has a zero at $z= -1$. The frequency response is $H(z)$ as z sweeps over the unit circle, thus giving a numerator magnitude as the difference between z at any point on the unit circle and the zero location: $z-q_z$ (or the multiplication of multiple such magnitudes if there is more than one zero), and a denominator magnitude given by the same for the pole locations: $z-q_p$. In this case the pole is at the origin, so $z-q_p=1$ for all $z=e^{j\omega}$. What also should be now clear is how the resulting phase response is formed since the net phase will be the difference between the phase of the numerator and the phase of the denominator (phases subtract in the division of complex numbers).

This type of nulling filter (zero-only) is not very effective given the gradual roll-off in frequency. To achieve very sharp nulls, place a pole very close to the zero; the closer the pole, the sharper the response! Given all poles must be inside the unit circle for a stable causal linear time invariant system, the magnitude of the pole would therefore be less than but close to 1.

This IIR approach is further detailed here: Transfer function of second order notch filter

Also this is an excellent write-up by Richard Lyons on linear phase nulling (or notch) filters that do provide a sharp notch with an FIR appraoch. This could similarly be translated to provide a notch at any frequency: https://www.dsprelated.com/showarticle/58.php


What is the application difference between extent and solidity in image processing?

I would suggest that Solidity is the better measure. Extent would be a cheap approximation. It’s cheap because computing the bounding box is cheaper than the convex hull (and really simple to implement).

Solidity is sometimes referred to as Convexity, but more commonly Convexity is the ratio of the perimeter of the convex hull to the perimeter of the object. These two measures are of course related, but not the same.

Solidity doesn’t quantify roundness. A rectangle has a Solidity of 1 (the maximum posible value). And so does a triangle. Roundness (or similarity to a circle) can be conveniently computed by comparing the area and the square of the perimeter, or the area to the product of the width and height (look for Feret diameters to compute these). Another approach is to determine the coefficient of variation of the distance of each boundary pixel to the geometric center of the object.

Solidity is useful to quantify the amount and size of concavities in an object boundary. Holes are also often included. For example, it distinguishes a star from a circle, but doesn’t distinguish a triangle from a circle.

Here is a good survey of object shape features: M. Yang, K. Kpalma and J. Ronsin, "A Survey of Shape Feature Extraction Techniques", in: Pattern Recognition Techniques, Technology and Applications, P.Y. Yin (Editor), I-Tech, 2008.

The DIPlib library implements the more meaningful of these features. (Disclosure: I implemented most of those features in DIPlib.)


Difference between 2$\pi f$ and $\omega$ in Fourier transform

You might be mixing two concepts, pertaining to (following Robert Bristow-Johnson)

the "analog context" (more formally, the "continuous-time case"),
the "digital context" (more formally the "discrete-time case").

The first concept corresponds to the continuous Fourier transform, for which you can use
a form of  normalized frequency cycles per second or Hertz ($e^{-j2\pi f}$), while the second is the angular frequency in radians per second ($e^{-j\omega}$).

The second one, with integer index $n$ relates first instance to discrete summations, e.g. the discrete-time Fourier transform (DTFT), where the Fourier kernel writes  $e^{-j2\pi f n}$ with $f\in \mathbb{R}$ (for instance) or $e^{-j\omega n}$ with $\omega\in \mathbb{R}$ (for instance).

The first concept is in use in the inverse DTFT, since it involves integrals.


What Is the Significance of a Large Residual When Applying Principal Component Analysis?

PCA dimensionality reduction seeks for a linear transformation, mapping the data to a lower dimensional space. The linearity comes from the mapping, which is simply a rotation (and translation if you take into account the de-meaning). In the end, it is nothing but a projection of the data.

So if your low dimensional embedding is not reconstructing the data truthfully, then this means that your data is not suitable for such linear mapping, i.e. not linearly separable. You should look into other techniques such as linear discriminant analysis (canonical variates), Kernel-PCA, which applies kernel trick to achieve non-linearity or manifold learning.


Bluestein's algorithm to evaluate the DFT from $f_o$ to $f_o + k\Delta_F$

The convolution in Bluestein's algorithm is a linear convolution, not circular. Since the approach is to use standard FFT's for the convolution operation, the underlying waveforms being convolved must be zero padded out to $M \ge 2N-1$ to compute the linear convolution accurately.


That's not an integral of variable $x$. The notation $x(a+T/2)$ stands for a function $x(\cdot)$ of variable $a$.
So applying the fundamental theorem of calculus, and assuming there exists a function $G(a)$ such that $G'(a) = x(a)$, then you will have :
$$ \int x(a+T/2) da = \int G'(a+T/2)da = G(a+T/2) + C $$
where the constant of integration, $C$, will be omitted in the definite integral :
$$ \int_{t+T}^{t-T/2} x(a+T/2) da = \int_{t+T}^{t-T/2} G'(a+T/2)da = G(a+T/2)|_{t+T}^{t-T/2}  $$
So the system has nothing with a square function.
Coming to its linearity, you can show this in line with the linearity of the integral operator...