Just despread the same way with two different integrate and dump times for your two symbol types with appropriate blanking on each between longer and shorter code sequence periods.

I would be tempted to do all the acquisition and tracking on the longer symbol and then simply despread the shorter symbols using the synchronization from the longer one— with that approach envision a receiver for the just the longer sequence with 62 chip gaps between each symbol.

Alternatively if each code repeats at the same location in each symbol, given the relatively short sequences, a sliding complex correlator using I and Q FIR filters with the sequence as the coefficients (in reverse order) using two samples per chip would be a very simple approach —- the complex I and Q output of the filter would have a spike in time after every alignment with those symbols (0 or 180 degrees once aligned depending on the data symbol and given I and Q would have phase to provide for carrier recovery).  Detect the spike with a threshold detector on the magnitude for acquisition, track the complex peak for carrier tracking and track the two samples on either side of the peak for early-late timing tracking. The peak sample will give you the same performance as an integrate and dump of the multiplication of the received sequence with the output of a code generator but given your peculiar waveform structure this may be a lot simpler.

Once all timing is established, a single 31 tap FIR filter on I can simply demodulate the other two symbols.

Summarizing this FIR approach:

One 126 tap FIR filter on I  for carrier and timing recovery/tracking and demod of 63 chip symbols.

One 126 tap FIR filter on Q for carrier and timing recovery/tracking.

One 31 tap FIR filter on I for demod of the 31 chip symbols.

As @WilliamHa points out in the comments, this approach will only work if the same code sequence appears in each symbol, starting at the same location each time (beyond being inverted to designate a "1" or "0" data bit), but will not work if the code generators run continuously causing the code starting location to shift by one sample at the start of each longer sequence, and at the start of the pair of shorter sequences, due to the relationship between the two codes.  In this case I would approach it using continuously running correlators where the outputs are blanked during the other symbols, using the longer code for tracking and carrier recovery and using the synchronized clock from that to demodulate the shorter codes. Similar to the FIR filter filter approach, there would be 6 correlators for the longer code: 2 for Early (I and Q), 2 for prompt (I and Q) and 2 for late (I and Q). Each of these should be running at at least 2 samples per chip. Only one I correlator is needed for the shorter codes at one sample per chip as the carrier and time reference is established from the longer code. The added complexity will be in timing the integrate and dump outputs of the correlator such that they are blanked during the shorter code sequences, but this would be determined in the timing recover process such that accumulation starts at the estimated start of the longer code and dumps after 63 chips, and then waits 62 chips before accumulating again. This structure assumes the code generators run continuously with the same chip rate, and the correlator structure is done by multiplying the code generator output with the received signal after having been down-converted to baseband. With this approach Carrier Recovery is done by monitoring the rotation of the I and Q prompt symbols of the longer code; a carrier error is determined using a complex conjugate multiply of one symbol with the previous and then multiplying that time I prompt to strip the +1/-1 data. Clock Recovery can be done by monitoring the magnitude of the Early and Late outputs of the longer code, where clock error = Early-Late. Acquistion is done by setting a threshold on the magnitude of the longer code.

Using the cross correlation property of the FFT, all the correlators could be done with FFT's as well, including one-shot joint frequency and delay acquisition if processing allows. I explain that approach further in this post: GPS signal acquisition


I expect that when this processed series are submited to denoising / thresholding will result in a clean series with smaller values than the submited signal

No, you can get values that are greater. For example, consider the Fourier series of a signal. The Fourier series basis functions are wavelets. If we approximate the signal with only a few of the Fourier series coefficients you will get values that are above the original signal.

E.g. from the (image from wikipedia)


Help filtering a logic level signal without degrading transitions

I've done a lot of work with things like this, hopefully the following will be of some help:

Since you've noted that the computation doesn't have to be performed in real time, this gives you plenty of options.  A natural choice, and I see someone has already suggested it, would be a median filter.  This usually works pretty well, so long as you choose the window length appropriately.  In your example that you showed, it looks like you may need to increase the window size (though it's hard to tell how many samples you have).  The main attraction to a median filter is that they will not smear edges, so long as you choose the window size appropriately; most standard smoothing methods will.

If you're over sampled, a legitimate option would be to just reduce your sampling rate and see where that gets you noise wise.

Other smoothing methods exist, but they'll typically smear your edges, and sort of ruin any edge/clock detection you'll want to do.  Since you have infinite computation time, median filters are probably your best bet in this realm.

Another option to check out would be a Savizky-Golay filter.  This is another data smoothing technique, but it attempts to model the signal using a polynomial.  MATLAB has a built-in that works great:

Savitzky-Golay Filtering in MATLAB

I've used this method on varying types of noisy signals, and if you're not constrained to running in real time, you can certainly tweak the parameters and hopefully find a combination that works for you.

Once you've got your data smoothed, a simple threshold would work I imagine.  You could statistically set this threshold based on the mean of the smoothed signal (this method assumes that the noise is fairly stationary).  Another method would be an edge detection filter, like say the first order derivative of a Gaussian or something.  You would convolve this filter with your smoothed signal, and then utilize some peak picking algorithm to get the start/stops of your transitions; this would be a little bit more robust than just utilizing some sort of differentiation filter due to the weights it implicitly applies.  You can get the filter in MATLAB by entering: diff(gausswin(sz,alpha)), naturally substituting the parameters for size and alpha that best work for you.

Another option that could forego smoothing with the right implementation is some type of hysteresis thresholding with a "hold" function built into it, given that you have some knowledge or expectation of what the values are when the the pulse isn't present.  This would involve having a double threshold, and creating some logic for deciding what values are above and what values are below.  The extra bit on top of a traditional hysteresis threshold would be that you'd need to specify some sort of holding mechanism, i.e. a segment is only valid if it breaks the top threshold for X many samples, and it is only invalid if it breaks the lower threshold for Y many samples.  This option would probably be the most work, as (from experience, I coded this for my own work fairly recently) it would take you a 100 or so lines of MATLAB code to implement.  Setting these thresholds is naturally the challenge here; I've found that a statistical approach is usually robust enough to work.

Last, but certainly not least, you could attempt to do some sort of adaptive filtering if you are able to characterize just the noise component of the signal, and said noise signal is stationary.  If you're able to capture just the noise in the signal, through an auxiliary sensor or some clever high pass filtering, you can use the Wiener-Hopf method to clean up the signal.  Here's an example:

By statistically characterizing the noise variance (not covariance in this case, since we only have one auxilliary signal) and cross-correlation, we're able to come up with an adaptive weight to remove the noise. For some background theory on this, consult Wiener Filter

Now at the signal is all cleaned up, a simple threshold would get you the edges.

If your noise isn't stationary, you can force it to be stationary in small windows through a process called time segmentation.  That's a bit involved though, and hopefully you won't need such a thing.

The MATLAB code that generated that plot is below:

up_s = zeros(1,250)+3;
d_s  = zeros(1,500)+3;
u_s  = ones(1,500)*10;
s = [up_s,repmat([d_s,u_s],[1,50])];
t = linspace(0,1,numel(s));
% set up the noisy signal
nz = sin(2*pi*500*(t));
snz = s+ nz;

% w = R^-1 r, from wiener-hopf method
% zero mean assumed for covariance and correlation;
% k terms cancel and are excluded
w = inv(nz*nz') * nz*snz';
% apply the weight to the signal
w_sig = snz - w'.*nz;

figure;
ax1 = subplot(2,1,1);
plot(snz);
title('Original Signal');
xlabel('Sample'); ylabel('Magnitude');
xlim([1 numel(snz)]);

ax2 = subplot(2,1,2);
plot(w_sig);
title('Signal Cleaned with Wiener-Hopf Filtering');
xlabel('Sample'); ylabel('Magnitude');
xlim([1 numel(snz)]);

linkaxes([ax1 ax2],'x');

tl;dr:
Median filtering is great, you just need to select the parameters appropriately for your sample rate. Savitzky-Golay filtering will also get you good results, given that you get the parameters right, which can be challenging. With any smoothing method, you can expect some smearing of the edges, but median filters really minimize this.

If you can characterize just the noise of the signal, you can use an adaptive filter technique like the Wiener-Hopf method; this method assumes that the noise is stationary, so that's something to keep in mind. With this method, you'll get essentially zero smearing of the clocking edges.


@ecook: Your time-domain notion of coherent integration is correct. (Where multiple periodic sequences must be time-synchronized before the averaging operation.)

The only way I could justify the remark that "coherent integration can be provided by DFT" is in the following scenario: let's say your DFT's input is an $x1(n)$ sinusoidal sequence and your DFT result is the complex valued $X1(m)$ sequence. Then some moments later you collect a new $x2(n)$ time sequence of the same sinusoidal signal and produce a new complex-valued DFT result of $X2(m)$. To average (integrate) those $X1(m)$ and $X2(m)$ sequences you would add their magnitudes. That's a form of incoherent averaging (incoherent integration), also called "scalar averaging." (Averaging $X1(m)$'s and $X2(m)$'s real and imaginary parts individually may lead to spectral energy cancellation and you don't want that to happen.)

If, on the other hand, your $x2(n)$ time sequence was time-synchronized with the previous $x1(n)$ time sequence you can add $X1(m)$'s and $X2(m)$'s real and imaginary parts individually. That's a form of coherent averaging (coherent integration), also called "time-synchronous averaging."

Know this however, as Peter Kootsookos points out, averaging phase angles is tricky. Have a look at Peter's explanation at:
http://www.dsprelated.com/showarticle/57.php


Signification of alternating signs of imaginary and real part of an FFT

This is exactly the expected result for a time delay: A delay in time is a linear phase in frequency; so you shifted the amount that would cause 180 degree phase shifts between your frequencies (Note that at your first frequency sample location you have a positive real and positive imaginary, which rotates to positive real and negative imaginary which is a -90° rotation, then follow the same process for your second frequency and you will see it is a +90° rotation; the phase change between the frequencies given is 180°). It is repeating at 180° due to the modulo property of phase, but you are really seeing the evidence of a linear increase in phase that is rotating past 180 degrees at each of your frequency samples where energy exists.

Consider a signal at given frequency going through a fixed length of cable; it will be delayed a certain amount in time, and given the frequency that time will relate to a certain phase shift between the signal at the input and output of that cable. Now consider a frequency that is doubled, the time delay is the same but there will be twice the phase shift between the signal at the input and output in this case.  This is the result of the Fourier Transform of a delay, the magnitude is one but the phase in linearly increasing at a negative rate.

So in your case we are seeing the resulting rotation of 180 degrees for the same distance between your frequencies tones (evidence of a linearly increasing phase).

If you plotted each frequency as a complex vector on an I-Q (real vs imaginary) diagram, as you change the delay of your signal in time, you will see the vector rotate accordingly with no change in magnitude. The higher frequencies will rotate proportionally faster.

Here is another view in case that was not clear describing the same effect and very specific to your case.  The top plot shows sine-waves in time for a first and third harmonic, and the delay at a fixed time offset such that the first harmonic has shifted -90° and the third harmonic as expected has shifted three times as much or -270°.

These last plots demonstrate your case exactly and shows the same thing on a complex plane which is consistent with your complex frequency outputs that have real and imaginary components. Note that I use the common notation I and Q for real and imaginary where "I" is "In Phase" and "Q" is "Quadrature Phase".

Here we see the same thing in that each frequency tone is a vector that rotates around the complex plane at a fixed rate; so the third harmonic is rotating three times as fast as the first harmonic (and the fifth is rotating five times as fast etc). After a given time delay, we see the diagram where the first harmonic has rotated 90° and the third has rotated 270°. To try to keep this clear, I first show the case where the first and third harmonic are in phase, starting at an angle of +45° (Such that I and Q are both positive as in your plots) and we see the exact condition where after a fixed delay the first harmonic has rotated -90° such that I is postive and Q is negative, while the third harmonic in that same time duration has I negative and Q positive. (-270° rotation.)

If we instead have a starting condition identical to your plots, where the first harmonic has both I and Q positive, but the third harmonic is 180° out of phase such that both I and Q are negative, we end up with the diagram below. We see the same thing where the first harmonic shifts -90° and the third shifts -270° after a given time delay. This is exactly the result from your plots: The first harmonic goes from I+, Q+ to I+, Q-, while the third harmonic goes from I-,Q- to I+,Q-.


how to detect and capture scene change in a recorded video

You may use Low Rank Models for Video Background Subtraction  as a pre processing.
Once you leave out the background even simple methods will do to discern any motion.


The first-order difference operation is a technique for numerical differentiation. It is the simplest method that I know of, and consists of just treating a discrete-time signal as piecewise linear between each pair of points and calculating the slope between them. The idea is that this time series of slopes should approximate the time derivative of the underlying signal.

$$
x_d[n] = \frac{x[n] - x[n-1]}{T}
$$

where $T$ is the sampling period of the signal $x[n]$.


How to Apply Gaussian Filter to Raw Audio Files in Python?

Do it using the following steps:

Look for the optimal frequency for the Low Pass Filter. Usually close as possible to the bandwidth of your desired signal.
Design an LPF filter according to the frequency above. If you're after Gaussian based filter you need to optimize the STD ($ \sigma $) parameter.
Apply the filter either using convolution, Using Numpy's convolve() function (Only in case of FIR Filter) or Scipy's lfilter() function (Which, in case of FIR Filter does convolution as well yet can also handle IIR Filters).

Remark
Gaussian based filters aren't optimal for the task you are after (Their passband isn't flat).
I'd would design a filter using Remez Method (See Scipy's remez()) which is Equiripple Method or using Chebyshev Filter Type II (cheby2()) and Butterworth Filter (butter()) which has flat passband (Namely they won't do any harm to the audio content).


You're describing DSL.

So, go and implement DSL! The common ADSL is an OFDM system.

  This being baseband is there a way of not having to use quadrature mixing to transmit the Real and Complex parts of the IFFT i.e. non-passband way of sending it?

Yes, you can just make sure not to use one half of your carriers. That's mathematically not different from using a half as long DFT and and shifting it up by half-bandwidth.

  I am trying to keep the demodulation part as simple as possible. Anyway I can get away with not using channel estimation i.e. non-coherent demodulation with DPSK or something similar?

I don't understand that question. OFDM is used exactly because you'd need an equalizer to use a broadband channel, and the channel estimator is OFDM, in essence.

Yes, there's OFDM systems that use a differential modulation on their subcarriers (namely, DAB/DAB+). The rationale, however, is usually a tradeoff between achievable rate and easy implementation: DAB has to deal with doppler spread, anyway, so it would require a significant amount of channel sensing if done non-differentially; if you're wasting much SNR on channel estimation, you can as well spend 3dB of BER curve and do differential modulation, especially if you're in the early 1990's and everything else would be computationally hard.

If you're doing fixed line communication, where the channel rarely changes, using a differential modulation is very inefficient compared to having a channel sensing preamble. You'd need to build a system with a higher bandwidth to achieve the same rate you could do with channel sensing, and that in a channel where you really want to avoid being forced to use higher frequencies.

  Any pointers to fully implemented systems that I can dissect to learn more?

ADSL.

  Considering the low data rate, I'll probably end up doing this on a low-end DSP.

... end up using existing DSL link technology instead of building this on your own.

I honestly don't think the modulation is the hard part here – the link and network layer becomes interesting. Twisted pair for a couple dozen meters sounds like you actually just want to use 10/100 Base-T. Modern network cards are surprisingly good at dealing with low-quality channels, and shielded twisted pair is actually a pretty nice channel (as opposed to what you claim, similarity to powerline).

If you really don't have the chance to implement half-duplex 10Base-T, you could probably simply go, physically, for simple serial RS-485; again, the problem is likely less a physical, but network layer design effort.


What is the largest "safe" order for the digital Butterworth filter of a given signal?

One cause is that higher order Butterworth filters have poles closer to the unit circle.  This nearby infinite gain point increases the likelihood of numerical instabilities.  (e.g. rounding/arithmetic/quantization noise may move a pole to the “wrong” side of the unit circle.)
Where the numerical noise will blow up depends on your executable code’s exact arithmetic sequence, and the quantization and number format used.  Often, dividing up a high-order filter into a sequence of biquad stages produces a more stable (better conditioned) sequence of arithmetic operations.


There seems to be a mistake in your derivations, the following code shows the correct one for $|\alpha| < 1$ ,  as indicated by @Matt.L $:

N = 128;         % simulation length (choose large enough)
n = -N : 2;      % simulation time span

C = 1.5;          % set some values
D = 2.78;
a = 0.33;         % a is for alpha

% x = -(n-2).*C.*(a.^(-n+1)).*stepfun(-n+2,0) + ...
%    -(n-1).*D.*(a.^(-n)).*stepfun(-n+1,0);

x = zeros(1,130);  % use the explicit code as stepfun() has problem with n.
for n = -N:1
    x(n+N+1) = -C*(n-2)*(a^(-n+1)) - D*(n-1)*(a^(-n));
end

M=2000;          % Select a FFT length for demonstration purposes
F=fft(x,M);      % compute the DFT of x[n] so obtained

w = linspace(0,2*pi-2*pi/M,M);              % choose DFT frequency samples
H = (C*exp(-1j*w)+D)./((1-a*exp(1j*w)).^2); % Evaluate the given DTFT

figure,plot(w,abs(H));title('H'); % compare the magnitudes of DTFTs
figure,plot(w,abs(F));title('F'); % Phases will be different as this
                                  % is an acausal sequence.

Derivation of the result for the left-sided sequence  $|\alpha|< 1$ case) is as follows:

Intuition plays a central for getting the answer as clean as possible which involves inverse Fourier transforms. In this example we shall use those three fundamental Fourier transform properties.

1- $x[n] \leftrightarrow X(\omega) \longrightarrow y[n]=x[-n] \leftrightarrow Y(\omega) = X(-\omega) $

2- $x[n] \leftrightarrow X(\omega) \longrightarrow y[n]=nx[n] \leftrightarrow Y(\omega) = j \frac{d X(\omega)}{d\omega} $

3- $x[n] \leftrightarrow X(\omega) \longrightarrow y[n]=x[n-d] \leftrightarrow Y(\omega) = e^{-j \omega d} X(\omega)$

Begin by $$x[n] = a^n u[n]  \longleftrightarrow X(\omega) = \frac{1}{1 - a e^{-j \omega} }$$

Let $y[n]=x[-n]=a^{-n} u[-n]$ then $Y(\omega)= \frac{1}{1-a e^{j \omega}}$

Apply second rule:
$z[n] = n y[n] = n a^{-n} u[-n]$ then $Z(\omega) = \frac{-a e^{j \omega}}{(1-a e^{j \omega})^2}$

Apply linearity:
$  w[n] =  -a^{-1} z[n] = - n a^{-1} a^{-n} u[-n]$ then $ W(\omega) = \frac{e^{j\omega}}{(1-a e^{j \omega})^2}$

Apply shift property on $w[n]$ to get rid of the numerator exponential.
$v[n] = w[n-1] = - (n-1) a^{-1} a^{-(n-1)} u[-(n-1)]$ then $ V(\omega) =\frac{1}{(1-a e^{j \omega})^2}$

Now apply shift property and linearity to $v[n]$ such that the resulting DTFT will be $$H(\omega) = \frac{ C \cdot e^{-j \omega} + D} {(1-a e^{j \omega})^2} = \frac{ C \cdot e^{-j \omega}} {(1-a e^{j \omega})^2} + \frac{ D} {(1-a e^{j \omega})^2}$$

which implies that
$$h[n] =  C v[n-1] + D v[n] $$ hence we get $h[n]$ as:

$$h[n] = C \left( - ((n-1)-1) a^{-1} a^{-((n-1)-1)} u[-((n-1)-1)] \right) - D \left( (n-1) a^{-1} a^{-(n-1)} u[-(n-1)] \right) $$

simplify to get the result:
$$ h[n] = C (-n+2) a^{-n+1} u[-n+2] + D (-n+1) a^{-n} u[-n+1] $$


How to perform a gaussian blur using FFT?

Here is the working code,

clear all; close all; clc;

I = double(imread('Cameraman.tif'))/255; % Divide only if image value in [0-255] range
figure,imshow(I);title('Cameraman original')
%I=rgb2gray(I);  % NOTE: This is only for RGB images, commented otherwise
sigma=3;
[M,N,s]=size(I);
f1=-fix(M/2):ceil(M/2)-1;
f2=-fix(N/2):ceil(N/2)-1;
[fx,fy]=meshgrid(f1,f2);

X=exp(-2*pi*pi*sigma*sigma*((fx/M).^2+(fy/N).^2));
If=fft2(I);
If=fftshift(If);
If=If.*X;
If=ifftshift(If);
I=real(ifft2(If));
figure,imshow(I);title(['filtered image with \sigma =', num2str(sigma)])


This is not the case. First of all, a second-order hold would use three sample points to compute an interpolation polynomial, but your suggested impulse response $\text{tri}(t)\star\text{tri}(t)$ is non-zero in an interval of size $4$ (assuming a sample interval of $T=1$, as you do in your question). However, the impulse response corresponding to a second-order hold must have a support of length $3$.

Now you could suggest that an $n^{th}$-order hold could have an impulse response which is the convolution of $n$ rectangular functions. In this case you would get the correct support size, but of course that's not sufficient.

An $n^{th}$-order hold computes a piece-wise interpolation using $n+1$ consecutive data points. This is in analogy with a zero-order hold using a single data point, and a first-order hold, which uses two data points. This definition is commonly used in the literature (see e.g. here and here).

It's straightforward to show that the second-order polynomial that interpolates three data points $y[-1]$, $y[0]$, and $y[1]$ is given by

$$P(t)=y[-1]\frac{t(t-1)}{2}+y[0](1-t^2)+y[1]\frac{t(t+1)}{2}\tag{1}$$

In order to find the impulse response achieving the interpolation given by $(1)$, we have to equate $(1)$ with the expression

$$y[-1]h(t+1)+y[0]h(t)+y[1]h(t-1)\tag{2}$$

If we choose the support of the impulse response $h(t)$ as the interval $[-1,2]$, which is equivalent to choosing the interpolation interval $[0,1]$, equating $(1)$ and $(2)$ results in the following impulse response of a second-order hold:

$$h(t)=\begin{cases}\frac12(t+1)(t+2),&\quad -1<t<0\\
1-t^2,&\quad0\le t \le 1\\
\frac12 (t-1)(t-2),&\quad 1<t<2\\
0,&\quad\text{otherwise}\end{cases}\tag{3}$$

The impulse response $(3)$ of a second-order hold looks like this:

I leave it up to you to show that this impulse response cannot be generated by convolving three rectangular functions with each other.


FIRFilter with low cut-off compared to sampling frequency

If the desired cutoff frequency is very low relative to the sample rate, then yes, you would need a very high-order FIR filter to implement it. However, in these sort of cases, it's usually a better idea to use a multistage approach.

For example, say you wanted to decimate a signal by a factor of 10000. If you implement this with a decimating lowpass filter, it would need to have a cutoff frequency less than $\omega = \frac{\pi}{10000}$. This would require a filter with tens to hundreds of thousands of taps, depending on how much attenuation you require in the stopband. This is a lot of processing, especially if you're trying to do real-time processing. What's more, with really long filters like this, you can run into numerical precision issues; the filter that you get might not even give you the frequency response that you designed it for!

Instead, consider an approach with two stages: two identical filters that each decimate by a factor of 100. In this case, your cutoff frequency must be less than $\omega = \frac{\pi}{100}$. And you'll only need hundreds to thousands of taps in order to implement each filter. The computational cost of two few-thousand-tap filters is much less than the cost of one filter with hundreds of thousands of taps. This only becomes more apparent when you notice that the first filter decimates by a factor of 100, making the cost of the second one negligible compared to the first stage.

In short, when large rate changes like this are desired, multi-stage approaches are almost always what you want. The same logic follows for interpolation as well. Lyons has an easy-to-understand discussion of this topic in his introductory text Understanding Digital Signal Processing.


The deconvolution operation in the code is just finding the impulse response of a filter made up of a numerator which is the signal to be deconvolved and a denominator which is (effectively) the filter to do the deconvolution.

As far as deconvolution algorithms go, it is a little simplistic.


Why averaging in Spectral magnitude domain not in complex domain to estimate spectrum of a process

You seem to understand the issues at hand.  In your situation, if the trial runs would be identical, except for the noise, averaging the signals (or the complex values) will definitely work, and I recommend that it how you do it.

The notion of averaging magnitudes (or magnitudes squared) across DFTs comes from when you are chunking a signal into multiple DFTs and want to aggregate the results.  In that case, if you average the raw DFTs you will get all sorts of cancellation/reinforcement effects depending on the frequency values of the signal compared to your DFT frame length.  When you average the magnitudes, the phase values are ignored and you will get more desirable results.

As is often the case wiht policy, it sometimes gets followed blindly without the understanding of the rationale behind the policy.  That may be what is behind some of the recommendations you have heard.

If you are trying to decompose your signal (assuming it is periodic) into its constituent harmonics, retaining the phase information is essential.

Hope this helps,

Ced


Do all real-time filters introduce phase lag, or can there be such a thing as "lag-less" or "zero lag" real time filters?

If a filter has any memory and that memory affects the current output, then that past data is lagged into the present.  If that filter is time invariant, then a memory of the present can potentially affect future output. e.g. is lagged into the future.  So, yes, unless the filter is memoryless or not time invariant.

Note that a unit impulse at the origin is memoryless, but impossible to implement in any finite size system due to the speed of light being finite.


Why is my resonance phase lag discontinous?

You seem to have 3 problems
1 - You have a phase wrap. It happens when the phase goes past 180 degrees. The phase will wrap back to -180 degrees. You can fix this by unwrapping the phase.
2 - Your phase is in normalized radians instead of being in radians or degrees. This is not a problem per se but when you try to unwrap the phase, you should be aware that phase unwrap blocks usually expect a phase in radians, not in normalized radians.
3 - The measured phase seems to be inverted compared to what you expect.


Your problem here is not a discriminative problem in which you have two (or more) classes you want to recognize. If it had been the case - for example discriminate infant cry with recording of bird songs, you would have compared each MFCC vector from the incoming audio stream to a massive set of frames either labelled as "bird song" or "infant cry" depending on the recording from which they came from. Note that the datasets would have been absolutely massive if your training sets had more than seconds of audio.

Your problem consists in measuring how similar your data is to a given category; or how well it fits a model trained on the example data - it is a "one-class" problem. In a Bayesian framework, this is equivalent to evaluating the probability of the "infant cry" hypothesis given the MFCC frame as evidence. I recommend you to frame your problem in Bayesian terms because it is likely that for your application, there are probably different costs involved in a false alarms or in missing the detection of a cry; so you'd like to work with a method which can balance them. Working in a probabilistic frameworks implies that you'll need to pick a model for the distribution of MFCC vectors on your training data. The simplest, but most computationally expensive approach is to use  Kernel density estimation - this is like a one-class, statistical cousin of k-NN (equally computationally intensive); but a more sensible approach given the amount of data would be gaussian mixture models.


What segmentation methods can be used for simple images?

One of the simplest method you can use is

"Region growing" scheme. Here the algorithm starts with a seed points, and keeps registering the pixel surrounded by it as long as the region as a whole has same property as before and after the inclusion of pixel. This is applied with reasonably good set of seed points recursively.

This property can be color variations, statistical properties such as mean and variance of intensity of the regions, or it can be texture. And the criteria for segmentation is homogeneity of the region.

The same principle but alternative implementation is region split and merge approach. Here a image is first divided into fixed parts - say 4 parts and corresponding scores of above properties is computed. Now the each sub region is split further, and we can check if any region is significantly different from others or same. If homongenety of sub parts is same essentially it is merged back. But if the element is realized to have different properties it now fall in different regions.

Read this for more details: http://vplab.iitm.ac.in/courses/CV_DIP/PDF/lect-Segmen.pdf


What you want to do is impossible with your technique. The phase of an FFT does not provide the details of an image. Even if it did, why would it enhance an image, if you exchanged the "details" of both images without any calculation?

The phase of an FFT of an image mainly holds the structure information, i.e. how the frequencies contribute to the image (in terms of phase offsets). This is valid for fine and coarse structures in the image. The magnitude tells you, how much each frequency contributes to the image. This information is in most cases less important, as can be seen by the general shape of the 2D magnitude spectrum (many strong low frequency contributions and a rapid decay towards higher frequencies).

I assume, what you would like do do is more a quite complex problem to be solved in image space. You might want to look for publications and algorithms on super-resolution for your application.


Recursively obtaining all-pole coefficients from complex cepstral coefficients?

If $h[n]$ is the impulse response of an all-pole minimum-phase filter of order $p$, then we have

$$\mathcal{Z}\{h[n]\}=H(z)=\frac{G}{1-\displaystyle\sum_{k=1}^pa_kz^{-k}}\tag{1}$$

where $G$ is some gain value, and $a_k$ are the predictor coefficients. Let $\hat{h}[n]$ be the complex cepstral coefficients of $h[n]$. If you know how to compute $\hat{h}[n]$ from $a_k$, then you're probably familiar with the following recursion:

$$\hat{h}[n]=a_n+\frac{1}{n}\sum_{k=1}^{n-1}ka_{n-k}\hat{h}[k],\quad 0<n\le p\tag{2}$$

From $(2)$ you can formulate a recursion for $a_k$:

$$a_n=\hat{h}[n]-\frac{1}{n}\sum_{k=1}^{n-1}ka_{n-k}\hat{h}[k],\quad 1<n\le p\tag{3}$$

with the initial value $a_1=\hat{h}[1]$.


Given Gain function, how to design a causal, stable and minimum phase IIR filter?

If $H(\omega)=e^{\alpha(\omega)+j\phi(\omega)}$ is a minimum phase frequency response, then the attenuation $\alpha(\omega)$ and the phase $\phi(\omega)$ are related by the following Hilbert transform relationship:

$$\phi(\omega)=-\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{\alpha(\Omega)}{\omega-\Omega}d\Omega$$

So $\phi(\omega)$ is uniquely determined by $\alpha(\omega)$. This is the corresponding wikipedia entry.


Question regarding AC power of ergodic process

Considering real valued WSS processes whose auto-correlation and auto-covariance functions are defined as
$$ r_x(\tau) = E[ x(t)x(t+\tau) ] $$

$$ c_x(\tau) = E[ (x(t)-\mu_x)(x(t+\tau)-\mu_x) ] = r_x(\tau) - \mu_x^2 $$

Then the following are just basic observations.

The total power of a WSS random process is given by
$$ P_x = E[ x^2(t) ] = E[ x(t)x(t) ] = r_x(0) $$
For a nonperiodic, WSS, RP to be ergodic in the mean, a necessary and sufficient condition is such that the auto covariance function $c_x(\tau)$ goes to zero as $\tau$ goes to $\infty$ :
$$ \lim_{\tau \to \infty}  c_x(\tau) = 0$$
Since auto correlation sequence and autocovariance sequences are related by
$$ r_x(\tau) = c_x(\tau) + \mu_x^2 $$

then the value of ACF as $\tau$ goes infinity yields the DC power as
$$ \lim_{\tau \to \infty} r_x(\tau) =\lim_{\tau \to \infty} c_x(\tau) + \mu_x^2 = 0 + \mu_x^2 $$

hence $r_x(\infty)$ is the DC power for an ergodic process.

The difference between the total power and the DC power is the AC power.

For a zero-mean RP, the total power and the AC power are the same.


adding phase noise in signal thorugh freuquency

First - this will depend on the nature of the phase noise you wish to model. If deltaf is a constant then it is just a frequency shift and really isn't considered phase noise. What you really need is deltaf to be a noise vector the same length as t. In matlab format you can do this by:

t=0:255;
fc=0.2; %just an example
deltaf = 0.001*randn(size(t));
B=cos(2*pi*(fc+deltaf).*t);

Notice that the frequency will have a finite variance in this case - but the phase variance will increase linearly over time - this may or may not be what you want. To have a finite phase variance you would do something like

t=0:255;
fc=0.2; %just an example
ph_n = 0.0001*randn(size(t));
B=cos(2*pi*fc*t+2*pi*ph_n);

Typically you can analyze these signals by using trigonometric identities and then assuming the noise portion is small or by using Taylor series expansions for the noise components.


Even though your calculation yields the correct result (in this case), the steps are not completely correct.

First of all, the covariance matrix of a random variable $n$ is given by

$$ C = E[nn^H]. $$

The nomenclature of variance is more suitable for scalar random variables, in my opinion (or it refers to the diagonal of the covariance matrix, i.e. the variance of each element of $n$).

Now, you state the covariance of $n$ should be $\sigma^2A$. Hence, we have

$$ E[nn^H] = \sigma^2A $$.

Now, you are interested in the covariance matrix of the random variable $y=A^{-1}n$, which is given by

$$
D = E[A^{-1}nn^HA^{-H}]
$$
Now, knowing that $A$ is constant, we can take it out of the expectation, to get

$$
D = A^{-1}E[nn^H]A^{-H}=A^{-1}\sigma^2A A^{-H}=\sigma^2A^{-H}.
$$

Now, knowing that $A$ is a covariance matrix and these are generally semi-positive definite and symmetric, we have $A^{-1}=A^{-H}$, so the final result is

$$
D = \sigma^2 A^{-1}
$$

which matches your guess. Note that in case you multiply $n$ with a different matrix $B$, the result is not $\sigma^2B^2A$ (as you might have guessed), but according to above calculation it is $\sigma^2BAB^H$.

I often refer to the following link, when it comes to expectations that involve matrices: http://www.psi.toronto.edu/matrix/expect.html


Curved noise floor and ugly unexplained peaks after FFT

Does the 'same configuration of your analyzer' mean the 'same sampling rate'? Is that equal to 41 MHz? It looks like you faced a well-known problem with a Nyquist criterion practical implementation. I mean that some time ago in digital oscilloscope ads there was an exact Nyquist ratio between a sampling rate and channel bandwidth (especially for cheap ones). Now it is about 10. (You may want to look at any ads of digital oscilloscopes) I remember there were several articles with figures similar to those posted here. To check this hypothesis just check a sampling rate of your spectrum analyzer (you wrote that artificial pikes are not seen in a spectrum analyzer with the same configuration.) or try to use a higher sampling rate in your data processing procedure if possible. In order to get at least an intuitive idea of the reason lies in the foundation of the problem perhaps you want to check this: Why is it a bad idea to filter by zeroing out FFT bins?
I suppose the following paragraph from the cited answer works for you: 'So if your original FFT input data is a window on any data that is somewhat non-periodic in that window (e.g. most non-synchronously sampled "real world" signals), then those particular artefacts will be produced' This phenomenon is called resolution bias error, or more commonly, the picket fence effect http://www.azimadli.com/vibman/thepicketfenceeffect.htm. Here you can find an illustration http://www.mechanicalvibration.com/More_on_spectral_leakge.html


How to transform a signal to go through specific points?

To have all the samples in the interpolated waveform land on the original samples implies using an integer interpolation rate (otherwise time domain distortion would need to be introduced, and I can't think of an application where that would make sense to do).
For integer interpolation, a very simple approach is to use the DFT with zero padding since the samples chosen will be exact, but this will have more deviation in between then other optimized interpolation methods further detailed below. This approach is very simple, so if the deviation between samples is of no concern then this approach would be favored and would be as follows:
Compute the FFT of the samples without zero padding.
Scale the result by $N$, where $N$ is the total number of original samples, zero pad the result and then compute the inverse FFT.  (Due to reciprocity the fft and ifft operations can be swapped).
In MATLAB or Octave this is:
N = length(test);
interp = 5;
result = ifft(N*fft(sequence), interp*N);

Other optimized interpolation methods using zero-insert and subsequent image filtering would also preserve the values, and with post processing the image filtering can be done with a zero-phase filter using filtfilt() and for an optimized solution in the least squares sense, together with multiband filter coefficients using firls(). The zero-insert and filtering approach is further detailed here: What is the impulse response used in an interpolation filter when upsampling?
If your samples are non-uniform to begin with, first resample to a uniform rate, where this post already addresses that approach: What is an algorithm to re-sample from a variable rate to a fixed rate?
But if the new samples do not fall on the same time locations (which will occur if there isn't an integer relationship in the duration between the original samples) then the interpolated samples should not be expected to match the original samples when they occur at new time samples that were not in the original sequence.


How can I determine the order of filter?

Note that the given slopes of $6n\textrm{ dB}/\textrm{octave}$ and $20n\textrm{ dB}/\textrm{dec}$, respectively, are asymptotic slopes for large frequencies. So you have to extrapolate the lines that are approximated by the magnitude responses at large frequencies.

For the given figure you could do the following. For the second order filter, look at the attenuations at $f=10^2\textrm{Hz}$ and at the next vertical line, which corresponds to $f=2\cdot 10^2\textrm{Hz}$. The ratio between those two frequencies is an octave, and the difference in attenuation is approximately $-22-(-35)\textrm{dB}=13\textrm{dB}$ (should be $12\textrm{dB}$).

For the $12^{th}$order system you could look at the attenuation values at $30\textrm{Hz}$ and $40\textrm{Hz}$, taking into account that the curve flattens towards its knee, and you need to extrapolate the straight line. With this in mind you would arrive at approximately $-31\textrm{dB}$ at $f=40\textrm{Hz}$, and at $-2\textrm{dB}$ at $f=30\textrm{Hz}$, which corresponds to a difference of $29\textrm{dB}$. Solving

$$\left(\frac{4}{3}\right)^x=2\tag{1}$$

for $x$ gives $x\approx 2.41$, i.e., you have to multiply $29\textrm{dB}$ with $2.41$ to obtain the attenuation per octave: $29\cdot 2.41=69.9\textrm{dB}$ (should be $72\textrm{dB}$).


$\sin(t)$ is no random process, because there's nothing random about it. You could add a random amplitude to get a random process:

$$x(t)=A\sin(t)\tag{1}$$

This is a random process because $A$ is a random variable. However, $x(t)$ is not stationary, but it is cyclostationary, i.e., its statistical properties vary periodically. You can make the process $x(t)$ stationary by adding a random phase:

$$\tilde{x}(t)=A\sin(t+\phi)\tag{2}$$

The phase $\phi\in [0,2\pi]$ is a uniformly distributed random variable that is independent of $A$. It can be shown that the statistical properties of $\tilde{x}(t)$ given by $(2)$ are independent of $t$, and hence, the process is stationary.


Frequency spectrum of an amplitude modulated signal

With $s$, $m$, and $c$ as you display, then the diagram of $S(\omega)$ is incorrect.  If it were, the spectrum of $m$ would be a triangle centered at zero frequency, which it is not.

$$
S(\omega) = \frac{A_cA_m}{4} \left[\\
\delta\left(\omega - (\omega_c + \omega_m)\right)\\
+ \delta\left(\omega - (\omega_c - \omega_m)\right)\\
+ \delta\left(\omega - (-\omega_c - \omega_m)\right)\\
+ \delta\left(\omega - (-\omega_c + \omega_m)\right)\\
\right]$$


signal reconstruction error in compressed sensing

Yes and No. If you consider your signal, noise free, meaning that all of the zero elements are really zero, then your answer is No, however though if the input noise is taken into the account, then the higher amplitude of the non-zero elements results in better signal reconstruction (high input SNR gives high output SNR).


How to reconcile "bandwidth" in the Shannon-Hartley Theorem with a spread-spectrum principle

If I'm only using bandwidth B1, doesn't that mean I can filter some of the noise out? And if so, would I be able to get S/N down? Wouldn't that be an alternative to occupying the full channel bandwidth?

Capacity increases linearly with bandwidth, but only logarithmically with SNR. So, increasing SNR by decreasing bandwidth is a bad idea if your noise isn't narrowband.

  Is spread-spectrum a bad idea to try to occupy the full channel bandwidth?

No, but there's more at play here:

Remember, $C$ is just an upper limit assuming you somehow (and it's not generally known how to do that) build the perfect transceiver including perfect channel coding.

Spread-Spectrum techniques can make for very robust transceivers, but they'd typically require $B_1 \ll B$. However, if $B$ is large, you'll likely be met with a non-flat channel, and that means things get a lot harder to receive correctly, and you'll end up building a complicated equalizer to counter ISI.

Especially fast-changing environments, which might be that the receiver, the transmitter, or potential paths change position or phase, getting a channel state estimate sufficient to do the equalization becomes hard enough to require you to add more pilots or otherwise redundant data to your transmission that the increase in bandwidth might not pay.
That's the main reason why for high-bandwidth systems (Wifi, cellular data), the world is leaving spread-spectrum behind in favour of multicarrier systems (typically, OFDM these days, but also FBMC): this has happened to WiFi (802.11b was DSSS, a/g/n/ac/p... use OFDM) and to Cellular standards (3G / UMTS was typically Spread Spectrum, 4G/LTE/4G+/5G is OFDM (and might turn out to do FBMC in the future for specific cases)).

Still, for smaller bandwidth systems, where the channel can be assumed to be tolerably flat, DSSS and alike play a very important role, because they allow you to send at a fixed PSD (say, $x$ dBm / MHz) over a larger bandwidth and thus allow for vastly improved SNR due to processing gain; in effect, you transmit the same bits, but you can put more power into them legally, because you're using a larger bandwidth. Yes, you infer the same increase in noise power; but noise doesn't correlate with your spreading pattern, so you'll have a net gain.

  it only allows you to increase SNR if you have narrowband noise.

Still haven't read that paper, so I don't know if you're misinterpreting the paper or if it's wrong (might simply make different assumptions):

Let's look at simple DSSS.

For every finite-power signal, we interpret its variance as power.

Let's assume we're dealing with a spread sequence $c_i$, $i=0,\ldots, L-1$, with $\left\lvert c_i\right\rvert = 1$. Let the symbols sent be $s$, $|s|=1$. This all works pretty well w.l.o.g., if you consider that the below also works with expected values instead of fixed amplitudes.

Then, what the transmitter sends for a single data symbol is a sequence $(sc_0, sc_1, \ldots, sc_{L-1})$. And at the correct sampling instant, the receiver calculates, from the receive samples $r_i$:

$$R = \sum_{i=0}^{L-1} c_i r_i \text.$$

In the noise-free case:

$$\begin{align}
R_\text{no noise} &= \sum_{i=0}^{L-1} c_i sc_i \\
&= s \sum_{i=0}^{L-1} c_i c_i \\
&= s L c_i^2 \\
&= s L 1^2 \\
&= s L
\text,\end{align}$$

which has the variance of $L^2$, since we scale $s$, which has a variance of $1$, with a constant $L$.

In the noise-only case:

$$\begin{align}
\text{Var}\left(R_\text{only noise}\right) &= \text{Var}\left(\sum_{i=0}^{L-1} c_i n_i \right)&\text{i.i.d}\\
&=L\text{Var}\left(n\right)
\text.\end{align}$$

So, the sum of noise has $L$ times the variance of the individual noise sample. Let's look at the SNR:

$$\text{SNR} = \frac{\text{Var}\left( R_\text{no noise}\right)}{\text{Var}\left( R_\text{only noise}\right)} = \frac{L^2}{L\cdot N} = L \text{SNR}_\text{unspread}\text.$$

So, I'd say, yes, spreading gain works clearly in your favor here.


I guess you are referring to an implementation like the one shown here. If you leave out the DFT/IDFT blocks, then you don't have the modulation step which is necessary for actually implementing a filter bank, and you will not be able to do appropriate sub-band processing between the analysis and synthesis stages.

When you test your analysis and synthesis stages without sub-band processing in between you will obviously get the same result if you leave out the DFT/IDFT blocks. Analysis and synthesis together are then a pure delay and you might as well replace everything by a pure delay, but this defeats the purpose of the filter bank which is to split the input signal into sub-bands for further processing, and then to reconstruct the original signal from these separate sub-bands.


Why does my digital upconversion result in two tones?

You can use the property of Fourier transforms which says that: $x(t)e^{j2\pi f_0t}\leftrightarrow X(f-f_0)$. That is, multiplying the time domain signal with a complex exponential will give you a shift in the frequency domain.

By multiplying by another sine, you get: $x(t)\text{sin}(2\pi f_0 t) \leftrightarrow X(f)*\text{FT}\big(\text{sin}(2\pi f_0 t)\big)=\frac{X(f-f_0)-X(f+f_0)}{2j}$, and so you see an extra peak at $120$ MHz.


Why LTI system cannot generate new frequencies?

One of the definitive features of LTI systems is that they cannot generate any new frequencies which are not already present in their inputs.
One way to see why this is so, comes by observing the output's Fourier transform $Y(\omega) = H(\omega)X(\omega)$
$$y(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau \longleftrightarrow Y(\omega)=X(\omega)H(\omega),$$ holds only when the impulse response $h(t)$  exists and it will exist only when the system is LTI.
From a little thought, guided by a simple graphical plot, and using the convolution-multiplication property above, one can see that the frequency region of support $R_y$ (set of frequencies for which $Y(\omega)$ is non-zero), of the output $Y(\omega)$ is given by the intersection of the regions of support $R_x$ and $R_h$ of the inputs $X(\omega)$ and frequency response $H(\omega)$ of the LTI system:
$$R_y = R_x \cap R_h$$
And from set algebra we know that if $A = B \cap C$ then $A \subset B$ and $A \subset C$ . That is, an intersection is always less or equivalent to what are being intersected. Therefore, the region of support for $Y(\omega)$ will be less than or at most equal to the support of $X(\omega)$. Hence no new frequencies will be observed at the output.
Since this property is a necessary condition for being an LTI system, any system that fails to posses it, therefore, cannot be LTI.


FIR design : compute coefficients for the same frequency response at another sampling frequency

There is no direct way of converting filter coefficients between two sample rates while maintaining the exact transfer function (at least where it's properly defined).
Re-sampling the impulse response can be done but will often result in extra latency, a longer filter, and some change in the frequency response.

In this case, however, you can derive the original filter specification from the filter coefficients and than re-design the filter using the same specification at a different sample rate. By visual inspection we can find that the ITU filter is an equiripple filter with

pass band edge: 20 kHz
stop band edge: 28 kHz
pass band ripple: 0.1dB
stop band attenuation: 40 dB

Let's say you want that filter at 4*44.1kHz instead of 4*48kHz. In Matlab you would do the following:

%% match the filter at 44.1 kHz
fs = 4*44100;
d = fdesign.lowpass('Fp,Fst,Ap,Ast',20000,28000,0.1,40,fs);
hd = design(d,'equiripple');
% convert to polyphase
h4 = 4*reshape(hd.Numerator',4,12)';

and the coefficients would be

    0.0057292     0.001552   -0.0015723   -0.0049864
   -0.0056096   -0.0017098    0.0050592    0.0099256
    0.0080942   -0.0011323    -0.012505    -0.017267
   -0.0093558    0.0088501     0.025828     0.027239
    0.0063161    -0.028386    -0.054115    -0.045155
     0.011432      0.10497      0.20287      0.26531
      0.26531      0.20287      0.10497     0.011432
    -0.045155    -0.054115    -0.028386    0.0063161
     0.027239     0.025828    0.0088501   -0.0093558
    -0.017267    -0.012505   -0.0011323    0.0080942
    0.0099256    0.0050592   -0.0017098   -0.0056096
   -0.0049864   -0.0015723     0.001552    0.0057292

However, there are a bunch of potential problems:

In this case we simply lucked out that the number of coefficients for the new sample rate was 48 as well. In other cases you may have to adjust the stop band frequency to dial in the desired number of coefficients.
This particular filter specification has been derived from the original sample rate, i.e. pass band and stop band are symmetrically spaced around the Nyquist frequency of 24 kHz. That is intentional. If your sample rate is different, you need to understand WHY it is different and whether the original specification is still appropriate. In the case of 44.1 kHz you may consider placing pass and stop band symmetrically around 22.050 kHz and maybe also consider lowering the pass band and/or increasing the number of coefficients.


What is the tightest known bound on the reconstruction error in compressed sensing?

According to Michael McCoy on LinkedIn:

The best bound for Gaussian matrices is in: W. Xu and B. Hassibi, Precise Stability Phase Transitions for l1 Minimization: A Unified Geometric Framework, IEEE Trans. Inf. Theory. Oct. 2011

See, for example, the curves in Fig. 1.


Why is modern control theory in the time domain instead of the S-domain?

I would say that the difference is not about time domain or frequency domain, is more about using transfer function or state space equations.
The state equation for a system with a single input and a single output is
$$\sum_{n=0}^{N} a_n \left(\frac{d}{d t}\right)^n y(t) = \sum_{n=0}^{N} b_n \left(\frac{d}{d t}\right)^n u(t)$$
The state space representation for a system
$$\dot{x}(t) = A \, x(t) + B \, u(t), \\ y(t) = C \, x(t) + D \, u(t)$$
where $y(t)$, $x(t)$ and $u(t)$ all can be vectors, extending this easily to multiple inputs and multiple outputs. For a high order system the number of parameters in state space representation is much larger, and there are multiple possible representations, however there are some canonical representations that are unique to a given system, and may have as few parameters as the transfer function.
Linking the two representations
This can be represented in the Laplace domain as
$$s \cdot X(s) = A \, X(s) + B \, U(s), \\ Y(s) = C \, X(t) + D \, U(s)$$
If you want to understand the dynamic of the system you can manipulate the state update equation as follows
$$s X(s) - A\, X(s) = B \, U(s) \\
 (s \mathbb{I} - A) X(s) = B \, U(s)$$
From the above equation you start to notice that when $X(s)$ is an eigenvector of You have $(s \mathbb{I} - A) v = s v - \lambda v$, and this will define one mode of the system $x(t) = x(0) e^{\lambda t}$, this is why you can study the dynamic of the system in terms of the eigenvalues of $A$.
Proceeding to the transfer function we can solve the above equation for $X(s)$
$X(s) = (s \mathbb{I} - A)^{-1} B \, U(s)$
and replace in the outptut equation
$$Y(s) = C \, X(s) + D \, U(s) \\
C \, (s \mathbb{I} - A)^{-1} B \, U(s) + D \, U(s) \\
(C \, (s \mathbb{I} - A)^{-1} B + D) U(s)
$$
And finally the transfer function is
$$H(s) = \frac{Y(s)}{U(s)} = (C \, (s \mathbb{I} - A)^{-1} B + D) U(s)$$
The result is so difficult to write that I will skip (and that is what is called classical control).
This equation will have rational functions (ratio of polynomials) in $s$, If a system has multiple inputs and multiple outputs it will be a matrix of scalar transfer functions mapping the effect of each input to each output.
Apart from notation compactness, linear algebra methods scale well to higher dimensions. One example of other advantage that I can cite is related to my master's degree. You can study system stability using Lyapunov stability condition, and this leads to linear matrix inequalities that can be applied to uncertain systems (with uncertainty in the model matrices).


The spread in the plot is the peak to peak phase, while the integrated phase noise is an rms quantity- The OP in his comment below this answer reported a measurement of 35° in the constellation and a 4.6° "integrated phase jitter". If this was rms phase then the measured peak as a ratio to the standard deviation would be $35\sigma/(2 x 4.6) = 3.8\sigma$ which is not an unreasonable result.

I had first assumed the OP was comparing an rms cycle-cycle jitter measurement in units of radians to the peak-peak phase deviation as shown on the constellation, hence the more detailed descriptions below that I will leave for future reference.

Determining RMS Phase

The rms phase error over the measurement duration can be determined from the SSB phase noise measurement $\mathscr{L}_\phi(f)$  using

$$\phi_{rms} = \sqrt{\int_{f=f_{low}}^{f_{BW}} S_\phi(f)df}$$

Where:

$S_\phi(f) = 2 \mathscr{L}\phi(f)$

$\mathscr{L}\phi(f)$: SSB Phase Noise (the values from the OP's plot above after converting from $\tt dBc/Hz$ to $\tt rad^2/Hz$)

$S_\phi(f)$: DSB Phase Noise (accounting for both upper and lower sidebands of phase modulation).

Th conversion from $\tt dBc/Hz$ to $\tt rad^2/Hz$ is $10^{k/10}$ where k is the value in $\tt dBc/Hz$.

With the range of integration further detailed below.

Determining RMS Jitter

A cycle-cycle jitter jitter measurement includes an effective high-pass filtering of the phase noise prior to integration, and converting to units of time.

The rms cycle to cylce jitter can be determined from the SSB phase noise measurement $\mathscr{L}_\phi(f)$  using:

$$t_{rms} = \frac{T}{\pi}  \sqrt{\int_{f=f_{low}}^{f_{BW}} (1-\cos(2\pi f)) \mathscr{L}_\phi(f)df}$$

This is arrived at by doing the following:

High pass filter $S_\phi(f)$ and integrate and square root to get rms phase:

$$\phi_{rms} = \sqrt{\int_{f=f_{low}}^{f_{BW}} |H(f)|^2S_\phi(f)df}$$

Convert units of phase in radians to units of time in seconds using:

$$t_{rms} = \frac{\phi_{rms}}{2\pi } T$$

Where

$H(f)= 1-e^{-j2\pi f T}$: Frequency response of cycle to cycle subtraction

$|H(f)|^2= (1-e^{-j2\pi f T})(1-e^{+j2\pi f T}) = 2-2\cos(2\pi f)$

$T$: Cycle time in seconds

$f_{low} = 1/T_m$ in Hz

$T_m$: total time duration of the data used for the rms jitter measurement in seconds.

$f_{BW}$: Channel bandwidth (all filtering prior to jitter measurement) in Hz.

This is further explained as follows:

The measurement of cycle to cycle jitter is done by computing the rms of the difference in time from the expected zero crossing location between adjacent symbols.

This measurement of cycle to cycle jitter is a high pass process - the jitter in one cycle is subtracted from the jitter in the previous cycle so is the filtering process given by the subtraction of the phase error (given in units of time) in one cycle from the phase error in the previous cycle.

As a discrete process this is $H(z)= 1-z^{-1}$ followed by an rms computation (compute the standard deviation of the time difference). Since the phase noise measurement represents a continuous time process, I used instead a continuous time integration with the equivalent filter $H(s) = 1 - e^{-sT}$. This is a completely valid approach since selecting the time error at each symbol location is a decimation process of the underlying continuous time function (of error vs time), which folds in all the upper frequency noise components, which then get captured in the integration.

This is actually a comb filter response, with a zero at DC forming a highpass over the measurement bandwidth given by the symbol rate.  Above cutoff and the noise is  uncorrelated so adds in power for a +3 dB gain, and below cutoff the correlation increases by the inverse of the frequency providing a $1/(f_c-f)$ rejection of the phase noise.  The easy way to see this is to consider the very very low frequency components of phase noise —- for those components the time shift in adjacent symbols is the same so is rejected when we compute the subtraction of the zero crossing offset between those two symbols. Below shows an example frequency response of the cycle to cycle measurement assuming a symbol rate of 10 MHz. Individual frequencies components above cut-off will have strong correlation every time the phase cycles through $2\pi$, hence all the nulls in the actual "comb" filter response.

Our observation time $T_m$ is given by the duration of the data set is another high pass filter (in order to see all the way down to DC you would need to observe forever!). An exact equation would add this additional high pass filter as an additional 20 dB/decade roll-off toward zero with a cutoff at $1/T_m$ and then extend the integration lower, but I have found it to not be necessary in all practical applications, getting the same result by simply starting the integral at $f_{low} = 1/T_m$ (similar to an equivalent noise bandwidth computation where the resultant filtering response is equivalent to a brickwall filter). So a brick-wall highpass with cutoff $1/T_m$ is used to model the observation duration, and similarly the channel bandwidth $f_{BW}$ is reasonably approximated as a brickwall lowpass filter since the phase noise is typically an insignificant contribution at the larger frequency offsets defined by the channel.

In summary, process the double-side-band (DSB) phase noise ($S_\phi(f)$) which is double the SSB phase noise ($\mathscr{L}\phi(f))$ with the high pass function that also has the additional noise power gain of +3 dB, and then integrate the result for total integrated power. Take the square root of this to get the standard deviation and this should match jitter once you convert between units of seconds as in jitter and units of radians as in phase noise using $t = \phi T/2\pi $ with $T$ as your symbol duration.


How do I analyze a music file ( only contains 5 seconds of note A) using the FFT plot to determine the note being played?

An FFT magnitude peak frequency (derived from a single FFT result index) is not reliably the same thing as the music pitch frequency or the musical note.  Note pitch is a psychoacoustic phenomena, and more reliably related to periodicity than the (sinusoidal) spectral peak frequency (e.g. a periodic/repeating waveform shape does not need to look like a sine wave).  See FFT Pitch Detection methods: Autocorrelation or other? for a list of some better methods for estimating note pitch frequency or the musical note.


Extracting common signal without knowledge about noise

Taking the average is a reasonable heuristic. The term “best” needs a lot more clarification because it needs a criteria.  Assuming that the signal is exactly the same on both channels is actually a strong assumption, so if you consider small perturbations to your model, say a time delay between the two channels, which is a lot more reasonable than perfect alignment, the character of the problem has changed.

There are techniques like ICA but they work well when they do and don’t when they don’t.

If your question is of a theoretical nature, there are answers to questions like is there a common signal in both channels using non parametric statistics.

If this is a practical question, then you should look at your data very closely and see whatever else you can reasonably assume.


Subcarrier assignment in MIMO-OFDMA system

For a given OFDM symbol, no subcarrier can be assigned to two users simultaneously. Each subcarrier is uniquely assigned to a unique user. Else you have to implement an Interference cancellation at receiver (which would require having channel knowledge of other users as well, meaning global CSI) or at the transmitter and extremely complex interference alignment/precoding has to be done.

In practical systems like LTE or 5G, thus one subcarrier in a given OFDM symbol is allocated to only one user for the reasons listed above.


Pinhole camera model from houdini parameters

This is a long topic to fully explain. I will try to write shortly, so please excuse the brevity.

Standard computer vision projection (ignoring distortion like Houdini) follows:

$$
\mathbf{x} = \lambda \mathbf{K}[\mathbf{R}\mathbf{X} +\mathbf{t} ]
$$

$\mathbf{R}$ is a $3x3$ orthogonal matrix, $\mathbf{t}$ is a $3x1$ translation vector. Camera position $\mathbf{C}$ is given by $\mathbf{C}=-\mathbf{R}\mathbf{t}$. $\lambda = x_z$ is a scale factor used for de-homogenization. The projection can also be compactified into a $4x3$ matrix:
$$
\begin{align}
\mathbf{x} &= \lambda\mathbf{P}[\mathbf{X}^T 1]^T\\
\mathbf{P} &= \mathbf{K}[\mathbf{R} | \mathbf{t}]
\end{align}
$$

Rotation is sometimes represented as a Rodrigues vector. To simply convert from one representation to another you might use cv::Rodrigues. Sometimes, (for some unknown reason!) people can also choose Euler angles for rotation representation. You should check which one Houdini uses, hopefully Rodrigues - if not you should compose $\mathbf{R}$ from 3 angles. To obtain $\mathbf{t}$ explicitly, just plug in $\mathbf{C}$ and $\mathbf{R}$ into $\mathbf{C}=-\mathbf{R}\mathbf{t}$ to solve for $\mathbf{t}$.

The only remaining unknown is then the $3x3$ camera instrinsic matrix, $\mathbf{K}$. In computer vision:

$$
\mathbf{K} =
\begin{bmatrix}
f_x & s & c_x \\
0 & f_y & c_y \\
0 & 0 & 1 \\
\end{bmatrix}
$$

The principal point $(c_x, c_y)$ is sometimes initialized to be $(\frac{w}{2},\frac{h}{2})$, which you can obtain from image size. Usually $s=0$. $f$ is just a bit more complicated:

Standard vision libraries (e.g. OpenCV) store the focal length in pixels, as in $\mathbf{K}$. In general, the relation to physical focal length is:
$$
f_{in-pixels} = \frac{w * f_{mm}}{w_{CCD-in-mm}}
$$
where $w_{CCD-in-mm}$ is the physical sensor size. In your case, I suspect that the aperture is $w_{CCD-in-mm}$, and focal length is the same as $f_{mm}$. If you plug these values you will get the focal length in pixels: $f=f_{in-pixels}$. Not sure about your units though.

Moreover, I assume that with $1$ x $aperture$ $= 41.4214$ you refer to a camera with square pixels with aspect ratio 1. Therefore it will not affect if you do the same calculation using height instead of width: $f_x=f_y=f$.

I hope this clarifies a bit. Let's keep refining this answer if you provide more info.


How to window a time domain signal for 6 symbols with cyclic prefix length of 512 samples?

Assuming that it would suffice to have an average result of all six symbols, I suggest that you simply concatenate all symbols with cyclic prefix in time and then window and FFT the whole dataset. This would give you the same result as an average of the result you would get for each symbol + prefix independently. Therefore the window length would be (2048+512)*6 = 15,360.

You would only consider overlap techniques when you need to continuously stream an fft process, or if you have a dataset that is not feasible to process with one fft. 15,360 is quite reasonable to do as one fft process.

Windowing involves multiplying the time domain data by the window weight sample for sample. After windowing the resulting windowed data is FFT'd. Purpose is to reduce spectral leakage from strong tones, thereby increasing the dynamic range and ability to see weaker signals together with stronger signals. This is at the expense of frequency resolution, or the ability to discern between two closely spaced tones that can be approximately the same power level (with decreased resolution, they would appear as one tone).

Note that the window will reduce the total power of the signal. For the measurement of power in bins where there are signals of similar power in adjacent bins (such as noise or waveforms that do not occupy single frequency bins), the non-coherent gain of the window would be used to compensate the measurement). For measurement of a single tone that occupies one bin (such as done with a calibration measurement), the coherent gain of the window would be used. The difference (in dB) between the coherent gain and non-coherent gain is $10Log10(RBW)$, where RBW is the resolution bandwidth. This makes sense as with a rectangular window (no further windowing other than selecting the data) the RBW is 1 bit, but any other window increases this bandwidth, thus if you sum square the bins to get the total power in the waveform, if you windowed other than the rectangular window you would be over-counting if there is similar power in adjacent bins since each bin includes the power in adjacent bins. The second link below goes into that further.

More on windowing:

How to apply Hamming Window?

How to calculate resolution of DFT with Hamming/Hann window?


Why does resampling change the frequency range?

Uniformly sampled pulse code modulation (PCM) audio such as WAV files only encodes frequencies up to the Nyquist frequency, which is half the sampling frequency. So for a sampling frequency of 500 kHz you can have frequencies up to 250 kHz, and for a sampling frequency of 192 kHz you can have frequencies up to 96 kHz. This is reflected in the range of the Y labels.


Why does White Noise in images imply noise in adjacent pixels are independent?

How did the professor state this? Did they say that the noise is iid?  That means that the samples are independent, identically distributed.  That means they are independent by definition (they are assumed to be so).

White noise, while a little more vague, also makes the assumption of independence. Again, this means that the noise is assumed to be independent.

However, if all that is said is that the samples are uncorrelated with the given covariance matrix, then you are correct: this does not imply independence.


Why does the Fourier Transform of the impulse look so different from the Fourier Transform of the impulse train?

This answer shows that $\sum\limits_{n=-\infty}^{\infty}
  e^{j  n \Omega T}$ is actually the Fourier Series of $\frac{2\pi}{T} \sum_{k=-\infty}^{\infty} \delta(\Omega - k\tfrac{2\pi}{T})$.

\begin{equation}
 \delta(t-nT) \frac{2\pi}{T} e^{j  n \Omega T}
\end{equation}
which means
\begin{equation}
 \sum\limits_{n = -\infty}^{\infty} \delta(t-nT) \longleftrightarrow \sum\limits_{n = -\infty}^{\infty} e^{j  n \Omega T}
\end{equation}
Now consider the exponential Fourier series of $\frac{2\pi}{ T} \sum_{k=-\infty}^{\infty} \delta(\Omega - k\Omega_s)$, because this function is periodic, so we can say
\begin{equation}
 \frac{2\pi}{ T} \sum_{k=-\infty}^{\infty} \delta(\Omega - k\Omega_s)
 =
 \sum\limits_{n=-\infty}^{\infty}
 c_n e^{j  n \Omega T}
\end{equation}
where
\begin{equation}
 c_n = \frac{1}{\Omega_s}\int\limits_{-\frac{\Omega_s}{2}}^{\frac{\Omega_s}{2}}
 \frac{2\pi}{ T} \sum_{k=-\infty}^{\infty} \delta(\Omega - k\Omega_s)
 e^{j  n \Omega T}
 d (\Omega T)
\end{equation}
and $\Omega_s = \frac{2\pi}{T}$ by definition.

Take $\alpha = \Omega T$,
\begin{equation}
 c_n = \frac{1}{\Omega_s}
 \frac{2\pi}{ T} \sum_{k=-\infty}^{\infty}
 \int\limits_{-\frac{T\Omega_s}{2}}^{\frac{T\Omega_s}{2}}\delta(\frac{\alpha
}{T} - k\Omega_s)
 e^{-j  n \alpha}
 d \alpha
\end{equation}
The term $\frac{1}{\Omega_s}
 \frac{2\pi}{T}  = 1$ and take now
\begin{equation}
 \beta = \frac{\alpha}{T}
\end{equation}
\begin{equation}
 c_n =
  \sum_{k=-\infty}^{\infty}
 \int\limits_{-\frac{\Omega_s}{2}}^{\frac{\Omega_s}{2}}\delta(\beta - k\Omega_s)
 e^{-j  n \beta T}
 d \beta
 \tag{1}
\end{equation}
We know that
\begin{equation}
 \int_{a}^b \delta(x - x_0) f(x) \ dx
 =
 f(x_0)
 \quad
 \text{if } a < x_0 < b
\end{equation}
and zero elsewhere. So the integral in $(1)$ contains one value within $[-\frac{\Omega_s}{2},\frac{\Omega_s}{2}]$ i.e. for $k = 0$, so
\begin{equation}
 c_n =
  \sum_{k=-\infty}^{\infty}
 \int\limits_{-\frac{\Omega_s}{2}}^{\frac{\Omega_s}{2}}\delta(\beta - k\Omega_s)
 e^{-j  n \beta T}
 d \beta
=
e^{-j  n (0) T}
=
1
\end{equation}
So $c_n = 1$ for all $n$. We have just proved that
\begin{equation}
 \frac{2\pi}{ T} \sum_{k=-\infty}^{\infty} \delta(\Omega - k\Omega_s)
 =
 \sum\limits_{n=-\infty}^{\infty}
  e^{j  n \Omega T}
\end{equation}
So whether you say $\sum\limits_{n=-\infty}^{\infty}
  e^{j  n \Omega T}$ or $\frac{2\pi}{ T} \sum_{k=-\infty}^{\infty} \delta(\Omega - k\Omega_s)$, it is the same thing.


A high-pass filtered signal has lost its low frequencies, so you can refocus the spectrum on a higher frequency range.
When you rectify it, the first phenomenon is that the signal becomes positive. Hence, it is not zero-mean anymore, and thus recovers low-frequencies that were attenuated in the high-pass filtered signal.

A second is the presence of higher frequency terms. Typically, when a sinusoidal component has period $T$, its absolute value doubles its period and generates higher-order harmonics.


The Hibert transform will produce a signal with 90 degree phase shift.  Look in Numpy documentation at scipy.signal.hilbert as a start.

I'm not sure what you are actually trying to do, but you mentioned 90 degree phase shift and radio in the same post, so you may find the following article of use:

http://home.comcast.net/~szemengtan/LinearSystems/hilbert.pdf


Incorporating Radial Distortion in Imaging Model

If you know the radial distortion parameters (e. g. by calibrating the cameras), then you should simply compensate for the distortion before you do structure from motion.  You can either undistort the images before you do anything else, or you can undistort the coordinates of the points that you are trying to match.


OpenCV: Problem at creating a Skeleton of Lines

You are actually drawing the skeleton of the background (brighter region).

change

cv::threshold(image, image, threshold, 255, cv::THRESH_BINARY);

in your codes to

cv::threshold(image, image, threshold, 255, cv::THRESH_BINARY_INV);

You should be able to get the right skeleton.


Camera calibration using checkerboard instead of 3D calibration object

Now with a checkboard given one image, I obtain the internal camera matrix and external kamera matrix.

No, you don't. You get a single coordinate transformation matrix that contains only the ratio between focal length and object distance. You cannot split that into an internal camera matrix (focal length) and external matrix (object distance) without further information.
If you want that extra information from only a single image, then the points of your calibration target must not be all in a single plane. This is why people use cuboid calibration targets. Alternatively, you can take at least two images of checkerboards from different perspectives, which also allows you split focal length and object distance.

Accordingly, am I right that it is necessary to do both, the calibration with the checkerboard to get an accurate estimate of K and then the calibration with the 3D calibration object to obtain the extrinsic camera matrix, so that the 3D calibration object defines the world coordinate system?

No, this is not necessary. A single image of a cuboid calibration target is enough to determine both the internal and external camera matrix and a single image of a checkerboard gets you nothing.


Verification of all steps needed to extract MFCCs from raw audio data

Point 7: Multiply and add the FFT with the MelFilter Banks. Therefore, if you are using 26 Filter banks , we will end up getting 26 dimensional vector at the end.
Point 9: Take the DCT of the 26 dimensional vector.
Point 10: Keep the first 12 coefficients and discard the rest.

I understood MFCC features by reading this link: http://practicalcryptography.com/miscellaneous/machine-learning/guide-mel-frequency-cepstral-coefficients-mfccs/

It is very simple to understand from the link above.


Why is my fundamental frequency at the 15th harmonic order

Yes you are correct that you are indeed seeing the equivalent of what would be the 15th harmonic of the fundamental frequency of the equivalent Fourier Series Expansion of the resulting equivalent continuous time waveform that was truncated in time to 40 samples (2 ms). According to the Fourier Series Expansion, $T$ would be 2 ms and the fundamental frequency for any continuous time domain trajectory over duration $T$ would be $1/T$. That fundamental frequency need not be the largest frequency nor need it be non-zero. The sum of that frequency with it's initial phase and amplitude with all the higher integer harmonics will converge to the original waveform.

That is a rather convoluted but mathematically accurate view. You are also seeing the expected result of the Discrete Fourier Transform for what must be a dominantly sinusoidal waveform that has a frequency of 15/40 of the sample rate (at or near 7.5 KHz). You are not sampling high enough to uniquely see higher order harmonics as they would be folded into the first Nyquist Zone from DC to $f_s/2$. (For example the 2nd harmonic at or near 15KHz would fold to 5KHz given your 20KHz sampling rate, or bin #10 in the plot above, but the signal may also have phase or amplitude variation causing such a sideband that is 2.5 KHz away- once sampled you can't distinguish these two effects).  You could certainly rotate (frequency translate) your signal to move to a lower bin as in your desired plot, but that will not resolve the ambiguity of the harmonics unless you do that frequency translation prior to sampling (as well as low pass filter out any signals that may still exist above $f_s/2$.)


What should be the singal duration in time doamin to correctly measure its frequency using FFT

For separating sinusoids (resolution of multiple peaks with clear gaps between them in an FFT magnitude result plot), you normally need an FFT length approximately double the reciprocal of the difference in frequency between the pair that is closest together.  This is roughly double the length over which 2 sinusoids are orthogonal, and is required in order to see the gap between, and to have 2 clear peaks instead of just one fat one in an FFT result plot.  This length increases if a non-rectangular window is used.


It's very hard to associate the numbers you're calculating with the actual signal power at the antenna. The SDR will filter, amplify, and AGC the signal, and then it'll convert it to digital. At best, you can say that the power you're calculating is somewhat related to the true signal's power.

The power you're calculating is in linear, not logarithmic, scale. To convert to decibels, use the usual formula. However, be aware that the reference inherent when using decibels is not well specified in this case.


Initial value Problem with difference equation

Like a linear constant coefficient ODE the general solution can be written as the sum of a particular solution and the general solution of the homogeneous equation $x(n+1)+x(n)=0$.

Because the right hand side is a polynomial we look for a polynomial particular solution, and a bit of experimenting will show this to be a quadratic:$$x_p(n)=n^2/2-n/2$$

Also the general solution to the homogeneous equation is obviously: $x_h(n)=k(-1)^n$, and the solution to your problem is the sum $x_p(n)+x_h(n)$ with $k$ chosen to satisfy the initial condition.


Minimum number of data points needed for a DFT to avoid spectrum leakage?

What you're actually asked to do is find the periods of the discrete-time sinusoids, and choose a DFT-length, such that all sinusoids occur an integer number of periods inside that window, hence no spectral leakage will occur.

I leave it as an exercise for you to determine the periods of the sampled sinusoids, but the end result of the DFT looks like this:

You can see that all frequencies are clearly resolved.

The problem (mentioned in the comments) of one sinusoid being located at exactly the Nyquist frequency (half the sampling rate) is not relevant for this toy example, even though in practice it is. What it means is that that sinusoid cannot generally be reconstructed uniquely. However, in your DFT you can see that there's a component at that frequency, and that's all they asked for.


does ifft of higher order equal upsampling signal

The longer FFT simply pads zeros to the original time domain signal. The inverse FFT should then be the original signal with zeros appended to the end.
The OP is correct that taking a longer FFT interpolates samples in frequency, and specifically provides more samples of the Discrete Time Fourier Transform. Why this occurs is clearer by observing the difference between the DFT and DTFT formulas as given below:
DFT:
$$F(k) = \sum_{n=0}^{N-1}x[n]e^{j k \omega_nn }$$
Where $\omega_n = 2\pi /N$
DTFT:
$$F(\omega) = \sum_{n=-\infty}^{\infty}x[n]e^{j \omega n }$$
For $\omega$ continuous from $0$ to $2\pi$ (for unique span)
Notice the two key differences: The time increment n in the DTFT extends to plus and minus infinity, and the frequency domain variable $\omega$ is a continuous function. As we add zeros to x[n] in the DFT, we approach the DTFT and ultimately for the same $x[n]$ all results of the DFT as more and more zeros are added will be samples on the continuous DTFT function.
In general related to this it is worthwhile to understand a universal property of all the Fourier Transforms related to signals that repeat or not (periodic or aperiodic) and continuous or sampled signals:
If a signal is continuous in one domain (time or frequency) it must be aperiodic in the other domain. And the counter-case is also true; if the non-zero values of a signal are discrete in one domain, it must be periodic in the other if we allow the domain to extend toward infinity. We see that immediately with the premise of the Fourier Series Expansion in that any single-valued analytic waveform over a finite interval can be represented by a discrete and infinite number of frequencies; adding together all those frequencies in the time domain over an infinite time duration will result in a periodic waveform.
With that observe the Continuous Time Fourier Transform (CTFT) as shown in the slide below. Regardless if the waveform itself is zero valued over a portion of the interval, the key point is the Transform itself uses a time axis that extends from minus infinity to plus infinity, and for that reason the frequency domain MUST be a continuous function: Frequency resolution is the inverse of observation time. If we have infinite time to observe a stationary signal, then we will have infinite frequency resolution. So the CTFT is of a signal in time that is continuous and aperiodic. Therefore the result in frequency is aperiodic and continuous.

Now compare this to the Fourier Series Expansion (FSE) as shown in the slide below, of a time slice of the same time domain signal. As previously explained, if the time is extended beyond the time interval from $0$ to $T$ the resulting waveform will be periodic. The FSE is specifically done over the interval from $0$ to $T$, but periodicity is implied as the exact same result in frequency would occur at the non-zero values if time was extended to minus and plus infinity. The frequency axis for the FSE is discrete since the results are defined for $n \omega_o$ with integer n. What also would occur if time was extended to infinity as explained for the CTFT is the frequency axis would become a continuous function, but the new interpolated values between the original samples would all be zero valued. The FSE of a signal in time from 0 to T is the same as a signal that is periodic in time over that same period. The signal in time is continuous and periodic, therefore the result in frequency is aperiodic and discrete (non-zero values are discrete).

Now compare this to the Discrete Time Fourier Transform (DTFT) as in the slide below for the same time domain waveform.  The only difference from CTFT is that the signal is now sampled in time. Thus the signal that is now discrete and aperiodic in time will be periodic and continuous in frequency: The same frequency result in the CTFT becomes periodic in the DTFT.

Finally the Discrete Fourier Transform (DFT) combines both effects of the FSE and DTFT; time is discrete as in the DTFT and restricted to a finite interval as in the FSE, with implied periodicity in that we can equally extend time to plus or minus infinity without changing the non-zero results in frequency. Thus the DFT for the same time domain waveform as the cases above would appear as in the slide below. The signal that is discrete and periodic in time would be periodic and discrete in frequency.


Drawbacks of upsampling using polynomial interpolation

Linear and more generally  polynomials are pretty common methods for interpolation or extrapolation, easy to implement, and simple in a causal setting (prediction). Parabolic interpolation (in the Fourier domain)  or Savitzky-Golay filtering are practically useful examples, especially for peak-like signals.

Polynomials however can be  poor at modeling finite support data, and are subject to instability, oscillations or overshoot: they don't have finite bandwidth in general, so they don't cope with standard "sampled signals" assumptions.

Addition (triggered by Fat32's answer): when you fit a polynomial under a least-norm, potentially with a penalty term, the regression may get rid of noises or outliers, so there is indeed a filtering effect). I have experienced that, for piece-wise continuous signals, using several weighted polynomial interpolation in parallel, and choosing the best one locally can be pretty efficient, even-though it does not lead itself to easy theoretical derivations (CHOPtrey).


Looking for a better algorithm for finding the slope of an input signal

I hope I understand what your'e asking correctly

To get the minimum slope of the image I would run an calculate the slope at each point (should be fast enough) and than take the minimal one.

For example :

x = 1:.05:3000;
y = 100+x + 300*sin(2*pi*x/1000);
origin = [0 0] ; % or origin = [ x(1) y(1) ];
slope = (y - origin(2)) ./ (x - origin(1));
[ minslope , minidx ] = min( slope ) ;
angle = atand( minslope );
figure;plot( x,y ,'-');
hold on;
plot([ origin(1) x(minidx)] , [origin(2) y(minidx)] ,'-r');

Will give you this slope:

Best of luck!,

Elad


Confused about Hamming window and its relation with frequency response specification

Windowing has nothing much to do with your filter design/response spec's, however, it has impact on your transition bandwidth that also depends on the window length.Your filters frequency response is mainly dependent on design parameters such as sampling rate, maximum frequency, transition bandwidth, cut-off frequency and impulse response length or kernel length.Now coming to windowing, these are some techniques basically used for obtaining smoother and ripple free impulse response h[n], so before your multiply your impulse response to the input/incoming samples, you must window your impulse response. For e.g. If you have a sinc impulse response, Ideally we know that its amplitudes run infinite on both negative and positive ends, hence in practice we just taper or cut off the ends of the sinc impulse in order to make it a finite length, therefore, this causes ripples in the impulse response due to the Gibbs Phenomenon (for more info on Gibbs Phenomenon, check online), So to eliminate these ripples and get a smoother impulse response, we use these windowing techniques.And each windows vary, for e.g. Blackman windows have better stopband attenuation than hamming window. Between why have you specified two stopband and passband frequencies?. Anyways,hope this gives you some insight to windowing!


Why is there a time-frequency precision tradeoff

why can I as a human hear lower frequencies sharply, but the FT cant?

No you can't. Human pitch discrimination at low frequencies is actually not particularly good. The "sweet spot" for human pitch perception is between 250 Hz and 4kHz and it gets rapidly worse above and below that. When tuning a bass instrument you typically try to reduce the fundamental and the lower harmonics as much as possible and emphasize the higher harmonics a lot since this makes it a lot easier to determine the correct pitch. A violin slightly out of tune makes everyone cringe: a double bass slightly out of tune is hardly noticeable.
So in this regard the human auditory system is no better than the Fourier Transform.


This depends a lot on how you implement it.

A single biquad takes about 10 arithmetic operations. (To be precise
a Transposed Form II takes 4-5 multiplies and 3 adds, depending on how the gain management is done).
Arithmetic operation translates into clock cycles of your processor. That depends a lot on the efficiency of your instruction set and how good yo are at writing code. On a good DSP you can do a single  biquad in 4 clock cycles (ca. 2 ops per cycle) and a general purpose ARM with a decent ALU this may by 15-20 cycles. On a crappy micro controller with no native hardware multiplier this can take up to 50-100 cycles.
Frame based vs sample based processing. Reading data from an ADC, moving it into the processor, and the moving it to the DAC creates a lot of overhead. That's why most implementations will process a frame of data and not one sample at a time. That "amortizes" the overhead across the entire frame. A good processor will also allow to do this with HW support, i.e. it will use DMA channels to move the data from/to DAC/ADC to internal memory without the CPU having to do anything.
The chips inside modern smart phones (even the cheap ones) are quite powerful: they are typically quad core ARMs with some HW accelerators (Neon, SIMD) for media processing. They all clock in the multiple GHz range and have good HW support for peripherals and graphic accelerators. My trusty old Samsung S8 has a 8 core 2 GHz Snapdragon: https://www.zdnet.com/article/qualcomm-snapdragon-835-what-does-the-kryo-280-adreno-540-spectra-180-x16-and-hexagon-682-mean-for/#:~:text=The%20Snapdragon%20835%20is%20comprised,result%20in%20optimized%20power%20savings.

Let's a run a simple example: A well coded biquad on a decent ARM core should take about 15 cycles. For 8 biquads running at 44.1 kHz on a stereo signal (two channels) that would be just about 10M cycles/second. That's about 1% of a single core at 1 GHz or 0.25% of the entire chip. On my Snapdragon 835 it would be less than 0.03%


Does fading memory mean impulse response with finite support?

What you want is theoretically impossible: if $h(t) \neq 0$ only on a finite interval, then it is not bandlimited. Fortunately, all real-world signals and systems have the property that $H(f) \rightarrow 0$ as $|f| \rightarrow \infty$. This allows you to approximate $h(t)$ very closely with only a finite number of samples. For engineering purposes, this is often more than enough.
Note that the property above means that $h(t)$ has finite energy and is absolutely integrable (it is in $L_1$). For the math to work out, it is often also required that $h(t)$ is in $L_2$, or that its square is absolutely integrable. These are the terms that are commonly used to describe these kinds of signals and systems.


Fourier transform of a signal and its autocovariance function

The interest in looking at the Fourier transform of the autocorrelation (which is usually equivalent to autocoveriance) comes from the Wiener-Khinchin theorem that states
Consider a wide-sense stationary signal $x(t)$. Then the power spectral density $S_x(f)$ is given by
$$S_x(f) = \int_{-\infty}^{\infty}R(\tau)e^{-j2{\pi}f{\tau}}d{\tau}$$
Where $R(\tau)$ is the autocorrelation function of $x(t)$.
There are times where the Fourier transform of a signal does not exist. However, this result still allows us to analyze the signal in the frequency domain and retrieve information by using the autocorrelation and not the signal itself.
This is also helpful when analyzing LTI systems. Given a system $h(t)$ with input $x(t)$ and output $y(t)$, we have a relationship between the power spectral densities given by
$$S_y(f) = S_x|H(f)|^2$$
This aids in estimating properties of the signals as well as the system itself, depending on what information is already known.
In the discrete domain, the DFT can always be calculated for a finite-length signal. The theorem still holds however, and despite being able to "take" the Fourier transform as in the DFT, errors may begin to appear when estimating parameters.
References by OP's request:
https://en.wikipedia.org/wiki/Wiener%E2%80%93Khinchin_theorem
https://mathworld.wolfram.com/Wiener-KhinchinTheorem.html
https://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-011-introduction-to-communication-control-and-signal-processing-spring-2010/readings/MIT6_011S10_chap10.pdf


External/Internal Description of a System

External. That's the whole point: The behavior of any LTI system can be completely characterized by it's transfer function, impulse response, differential equation, etc. (which can all be converted into each other).

That's a "black box" description and gives no information on the internal details of the system. There many different physical realizations of the same system, but you can't tell from just the outside behavior


relationship between energy and sampling rate

You must observe the nyquist frequency when sampling a signal.
In order to sample a signal without introducing artifacts, you must first filter out everything in the signal that has a frequency higher than half of your sampling rate.  If you sample a 1000Hz, then you must first filter out everything above 500Hz.

Yes, changing the sampling rate can change the bandwidth of digital representation of the signal.  A lower sampling rate means a smaller bandwidth.

This will only improve crosscorrelation if the noise is all high frequency stuff, and the signal you are trying to detect is in the lower frequencies.  You could get the same effect, however, by using a low pass filter  on the sampled data.  That would also remove the high frequencies and make the correlation clearer.


difference between Gabor and log-Gabor function

The 1D gabor filter has the following form in the frequency domain:

$$G_{b(\sigma,\omega_0)}(\omega) = \text{exp}\left(-\frac{\sigma^2}{2}(\omega - \omega_0)^2\right)$$

The 1D log-gabor filter is:

$$G_{l(\sigma,\omega_0)}(\omega) = \text{exp}\left(-\frac{\ln^2(\omega/\omega_0)}{2\ln^2(\sigma)}\right)$$

Log-gabor filters are used because they have 0 DC component for arbitrary large bandwidth, and size distribution of features in an image is often logarithmic.

The excellent paper On the Choice of Band-Pass Quadrature Filters by Djamal Boukerroui, J. Alison Noble and Michael Brady explains further differences for these two filters as well as others.

Here is some MATLAB code for the 2D version:

[rows, cols] = size(I);

% FFT mesh
[ux, uy] = meshgrid(([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
    ([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2)));
ux = ifftshift(ux);   % Quadrant shift to put 0 frequency at the corners
uy = ifftshift(uy);

% Convert to polar coordinates
th = atan2(uy,ux);
r = sqrt(ux.^2 + uy.^2);

% Create spectrum
filterFFT = 1.0/wavelength;
filterFFT = exp((-(log(r/filterFFT)).^2) / (2 * log(sigma)^2));
filterFFT(1,1) = 0;

% Filter image
I_filtered = real(ifft2(fft2(I) .* filterFFT));

Advantage of 0 DC

The advantage of having 0 DC component is that the response of the filter doesn't depend on the mean value of the signal. It also means you don't have an infinite impulse response. The Gabor filter has a small non-zero component around DC that is usually removed.

For feature detection, say you designed a filter kernel to match a step (1D) or edge (2D). If the filter had a DC component then the response would vary with the mean image level. e.g. consider the sobel edge detector, and the same kernel with increased values.

I = double(imread('the_image.jpg'));
I = I(:,:,1);
f1 = [[-1 0 1];...
      [-2 0 2];...
      [-1 0 1]];
f2 = [[0 3 4];...
      [1 3 5];...
      [0 3 4]];
I1 = conv2(I,f1,'same'); images(I1); colorbar; pause;
I2 = conv2(I,f2,'same'); images(I2); colorbar; pause;
I3 = conv2(I+100,f1,'same'); images(I3); colorbar; pause;
I4 = conv2(I+100,f2,'same'); images(I4); colorbar; pause;

In this code f2 = f1 + 3. Only convolving with f1, which has 0 DC, gives the same result when the mean value of I is changed. The image responses 'look' the same, but only the 0 DC one is useful for feature detection.


If the transfer function between the mains signal and your contaminated signal is fixed forever, then you could cancel the interference with a sine wave that is locked to the mains, and adjusted in both amplitude and phase to cancel your interference. However, using an fft is not the optimal way to do this, as the true mains frequency will likely fall in between the bins, so you will have a frequency error. A better way is to use a dsp phase-locked loop, the details of which are too complex to describe here.
However, it is unlikely that the transfer function is fixed. Therefore you will need an adaptive algorithm to adjust the amplitude and phase of your cancellation signal. This adaptive algorithm will look at (input signal - cancellation signal)  and attempt to minimize the portion that is correlated to the mains frequency. Unfortunately the adaptive algorithm will be affected if there is a component in your desired signal that is close to the mains frequency, which will reduce the attenuation of the mains frequency and also add attenuation to the desired signal. So it behaves very much like a notch filter, so you might as well just use a notch filter as RBJ has recommended.
You can minimize this interference effect by making the adaptive algorithm very slow, but again this is roughly equivalent to using a higher-q notch filter, and both systems will take a long time to settle out.


First if y[n] and x[n] were periodic, then wouldn't the ratio of their periods always be rational since each period is always integer?

No. I mean. Yes.

$x(t)$ might be e.g. $2\pi T$-periodic, and $y(t)$ could be $5T$-periodic ($T$ being the sample period). The choice of using capital $N$ for the period was a confusing one, but signals don't need to have integer sample-multiples as periods.

BUT, and that's where I'm with you, if you say $x[n]$ is periodic, it means the sampled signal $x[n]$ is periodic – which it can't be, if the original sample would have a rational-amount-of-samples period.

And thus, you're right: if the sampled signals are periodic, they need to have a rational ratio in periods, and thus, their sum will always also be periodic.


A typical approach would be the Pan-tompkins algorithm. You can also apply a first order low pass filter or a mean filter to improve results.