How to solve ADMM for TV Minimization Problem For Different Sizes $A$ and $x$ in $Ax=b$

The Error in the Model

The problem is in the dimensions of the Linear Operator $ A $ in your model compared to the Data Matrix $ X $. The number of columns of the matrix $ A $ must match the number of pixels in each column of $ X $ (Each image). While in your case it matches the number of images.

Matrix Form of 2D Linear Operator

Let's try to build the model correctly. Assume our data (Images) is given by the set of 2D matrices $ {\left\{ {X}_{i} \right\}}_{i = 1}^{m} $ where $ {X}_{i} \in \mathbb{R}^{p \times q} $. Given a Linear Operator $ \mathcal{A}: \mathbb{R}^{p \times q} \to \mathbb{R}^{k \times l} $ and we have $ {B}_{i} = \mathcal{A} \left( {X}_{i} \right) $.

In order to create the matrix form we should represent the Linear Operator by the matrix $ A \in \mathbb{R}^{\left( k l \right) \times \left( p q \right)} $ which is the matrix which applies the operator on the Column Stacked (See Vectorization Operator) vectors. So we have $ \boldsymbol{x}_{i} = \operatorname{Vec} \left( {X}_{i} \right) $ and $ \boldsymbol{b}_{i} = \operatorname{Vec} \left( {B}_{i} \right) $.

By defining $ X = \left[ \boldsymbol{x}_{1}, \boldsymbol{x}_{2}, \ldots \boldsymbol{x}_{m} \right] $ and $ B = \left[ \boldsymbol{b}_{1}, \boldsymbol{b}_{2}, \ldots \boldsymbol{b}_{m} \right] $ we have $ A X = B $.

With this definition everything will work as expected. Same logic will work for $ F $.

Remark

If you want $ F $ to work along the rows of $ X $ then set $ Z = F {X}^{T} $.

ADMM for the Vector Case

First, pay attention that the ADMM works on vectors (You can work with Matrices but then you need to update the Prox operations accordingly.

So $ \boldsymbol{x} \in \mathbb{R}^{n}, A \in \mathbb{R}^{m \times n}, \boldsymbol{b} \in \mathbb{R}^{m} $ and $ F \in \mathbb{R}^{o \times n} $ where $ o $ is set to match the operation of TV on the image columns and rows which $ \boldsymbol{x} $ is a column stack of.

So now the terms do align propely: $ {\left( {A}^{T} A + p {F}^{T} F \right)}^{-1} \in \mathbb{R}^{n \times n} $ and $ \left( {A}^{T} \boldsymbol{b} + p {F}^{T} \left( \boldsymbol{x} - \boldsymbol{u} \right) \right) \in \mathbb{R}^{n} $. So the calculation of $ \boldsymbol{x}^{k} $ is well defined.

Given that $ \boldsymbol{z} \in \mathbb{R}^{o} $ and $ \boldsymbol{u} \in \mathbb{R}^{o} $ their calculation is also well defined.

I guess something in your code doesn't match. But if the code have the terms $ A \boldsymbol{x} - \boldsymbol{b} $ and $ F \boldsymbol{x} $ well defined then everything else will work.


Can you represent an audio signal in any other domains besides time and frequency?

Basically it's all about breaking down the information into several 'bits'.  The actual audio signal is a time varying 'value', however often it is useful to consider it in a different form.  As an analogy, consider the number 256:  depending on what you're doing with your numbers, it might be useful to treat the number as 200 + 50 + 6, or 16 + 240, or 16*16, or maybe as 2^8; there are an infinite number of ways to treat the number, and which one you use depends on what you're trying to achieve.

The 'frequency domain representation' is an example of doing the above break-down, but with a signal rather than a number.  In this case, you're representing the original signal as a sum of sinusoids, all with different frequencies, amplitudes, and phases.  If you add them all together, you get your original signal back.  You could alternatively choose to represent it in a different way, such as with wavelets, or any other approach that may or may not even have a name (yet), if it is useful for what you're doing to it.  Maybe you could break it into 3 second chunks, then rearrange the signal in each of those chunks to be monotonically increasing, and remember how you reordered them.  This sounds a bit ridiculous, but there is an approach to manipulating MRI images with this sort of reordering (it doesn't use the chunking part, but that's just an example of how you can do whatever you like).

One advantage of sinusoids is that, like Nathan Day says, it relates to the way our ears interpret the pitch of sounds.  But a more important reason is that sinusoids are complex exponentials, which are eigenfunctions of linear systems; that is to say that linear systems are much much simpler to analyse if you consider the inputs and outputs as sums of sinusoids.  That is the main reason why Fourier analysis is so widespread and important.

Short answer to your question: there are an infinite number of domains you can represent an audio signal in.  For another popular one, see wavelets.


Morphology - Why does the concept of opening and closing images exist when we already have erosion and dilation?

As you said closing is completely defined by a dilatation followed by an erosion with the same structuring element (SE) (respectively opening is an erosion followed by a dilation with the same structuring element).

The key point is that is must be the same structuring element and this is much more specific than an erosion followed by a dilatation with different random SE.

Closing and opening operations are defined just for convenience, because giving an image and a SE the result can easily be predicted and interpreted by the user.

  The result of this particular operation is similar to dilation [...]

Not really, dilatation is more like a thickening operation while closing fills holes and smooth borders accordingly to the SE (opening with a disk rounds angular borders for instance).

Try erosion, dilatation, opening and closing on simple shapes (square, disk, triangle, hollow shapes, set of close shapes) with simple SE (share, disk, lines) to better visualize the effect of each of them. you can also try erosion followed by dilatation with different SE to see the difference with an opening.

To sum up opening and closing are defined for convenience only, as it is a specific case of erosion and dilatation composition that have noteworthy properties.


What is the story behind the story about SIFT descriptor?

The descriptor obtained from a $64\times 64$ neighborhood of interest point at the obtained scale.

It will divide this $64\times 64$ region to $16\times 16$ patches which lead to 16 patches.

For each patch we calculate the gradients and then find the dominant direction of gradients(which has some details), then taking the dominant direction as the reference direction we will divide the 360 degrees to 8 angular region each has 45 degrees, then sum over the magnitude of each gradients which lie in each angular region.

We could consider this as distribution or 8 bin histogram of gradient direction (considering strong gradients has more information we have to use them with higher weight in calculation of distribution so we use their magnitude as their weight which leads to sum over their magnitude). Then we will normalize these histograms.

At the end for each patch we have an 8 bin histogram and we have 16 patches which leads to 128 number descriptor.

By finding dominant direction our descriptor also becomes rotation invariant.
By using gradients our descriptor becomes invariant regarding to baseline illumination and by normalizing obtained histograms our descriptor becomes invariant to the contrast of image.


How can I find this all pass filter's phase?

phase of a filter is the angle (or "argument") of the complex number that is $H(e^{j\omega})$ at each value of $\omega$ .  you may need to "unwrap" the phase by adding or subtracting integer multiples of $2 \pi$ .

BTW, to post math to these pages, surround the $\LaTeX$ expression with single dollar signs for inline math and double dollar signs for math expressions on their own line.  if you need a quick guide for LaTeX, my suggestion is the one that lives at Wikipedia: https://en.wikipedia.org/wiki/Wikipedia:Math .  i'm still learning how to do other markup here.


Perfect reconstruction in two-channel filter bank

$G_0(e^{j\omega})$ and $H_0(e^{j\omega})$ are ideal low pass filters with passband $[0,\pi/2]$ and stopband $[\pi/2,\pi]$. $G_1(e^{j\omega})$ and $H_1(e^{j\omega})$ are ideal high pass filters that are complementary to the lowpass filters, i.e., their passbands coincide with the lowpass filters' stopbands and vice versa.
Note that shifting the frequency axis by $\pi$ transforms a lowpass filter to a highpass filter and vice versa. I.e., $H_0(e^{j(\omega-\pi)})$ is a highpass filter, and $H_1(e^{j(\omega-\pi)})$ is a lowpass filter. Consequently, $G_0(e^{j\omega})H_0(e^{j(\omega-\pi)})=0$ and $G_1(e^{j\omega})H_1(e^{j(\omega-\pi)})=0$. So the second term $(111b)$ corresponding to the aliased input spectrum vanishes. We're left with the first term, which is constant because we add an ideal lowpass spectrum and its complementary highpass spectrum.


I would implement it differently altogether.
Since applying the Wiener filter is pretty "cheap" I would create an Image called mWienerFilter.
Then:

vNoisePixels = (abs(I -  mMeanImage) > noiseThr);
I(vNoisePixels) = mWienerFilter(vNoisePixels);

That would be the best thing to do in my opinion.


Can $\delta(t+\infty)$ be a legitimate signal?

It's possible, in mathematics, to complete real numbers with "infinite" values, with sound topological properties; for instance non-standard analysis or the Extended real number line (discussion at Math StackExchange).

However, in this context, for any standard real number $t$, the rule is to set $t\pm \infty = \pm \infty$ (see Arithmetic operations). So, using your notations, for every $t$, you would have [CAUTION ADVISED] $\delta(t+\infty)=\delta(\infty)$, to which I am not able to give a meaning, else than something being a constant evererywhere equal to $\delta(\infty)$. Of course, $\delta$ is a distribution, and should not be treated as a function, but defined as a functional operator, or via limits of functions.

However, so far, I have never seen

a (serious and useful) use of the extended real system in signal processing,
a reference trying to define such a Dirac at infinity.

For the latter, I fear (just intuition) that constructing this with test functions can be troublesome, since two different limits are required (location and amplitude). This is just an analogy, but $\delta(0)\times \delta(0)$ is, as far as I know, not defined.


matlab How to design lpf\bpf\hpf without builtin functions

You seem to be confused with the time and frequency domain. Your code for an LP filter is the inverse Fourier transform of an Ideal low-pass filter. So by getting a frequency response for the function you wrote you'll get an ideal LP response you want.

However, there are two mistakes:

1) Your frequency for the ideal LP filter needs to be $[-\frac{m-1}{2} : \frac{m-1}{2}]$

2) You need to change the value for the 0-frequency only if M is an odd number.

Now, you can use the function freqz to get your freqency response, so your final code would be like this:

%[y,Fs]=wavread('sound.wav');
%m = length(y); %30586
m = 30586;
Fc = 500;               %% frequency cut
n = -(m-1)/2 : (m-1)/2;            %% kernel length / samples
%Fs = Fs;    %22000      %% sample frequency / (nyquist - %2*m)
Fs = 22000;
Tc = 2*pi*(Fc/Fs);      %% theta c
Hm = sin(Tc*n) ./ (pi*n);
if (rem(m, 2) == 1)
    index = (m+1)/2;
    Hm(index) = Tc/pi;
end

[H w] = freqz(Hm, 1, 1024); % 1024 - or how many points you want for your Fourier Transform

% Now your H is your frequency response, and w is digital frequency in [0,pi].
% You can now plot the amplitude response:
figure,plot(w, abs(H));

And you have your Low-Pass filter. However, I'd like to suggest using the many ways to create filters in Matlab, they were made for a reason, and allow for much easier filter design :)

One more remark, I have to suggest that you know the difference between analog and digital frequency, as they have different notation. Analog frequencies you write in uppercase, e.g. $F, \Omega$, and digital frequencies in lower-case, $f, \omega$. Here's a version of the code which is a little prettier:

M = 30586;
Fc = 500;
Fs = 22000;
wc = 2*pi*Fc/Fs;

n = -(M-1)/2:(M-1)/2;
h = (sin(n*wc))./(n*pi);
if ( M-2*fix(M/2) ) > 0
        index = (M+1)/2;
        h(index) = wc/pi;
        display(index);
end

[H w] = freqz(h, 1, 1024);
plot(w, abs(H));

Cheers :)

EDIT:

When you use freqz, the digital frequency w can be somewhat hard to control, because it's not created in the same way, or more clearly, the points are not selected like a normal frequency we create with, for example, F = 0 : Fs/(M-1) : Fs, even though it should be the "same" frequency.
The result is that it's not really possible to get the output filtered signal with Y = H .* X, and  y = ifft(Y).

However, I would suggest using the function filter(b, a, x) for filtering your signal, which works good with freqz. You can proceed to do just that:

% ------
% Create a random signal x
n = 0:M-1;
F = -Fs/2:Fs/(M-1):Fs/2;
Fx = [2000 7000 10000];
x = cos(2*pi*Fx(1)/Fs*n) + cos(2*pi*Fx(2)/Fs*n) + cos(2*pi*Fx(3)/Fs*n);
X = abs(fftshift(fft(x, Nfft))); % fftshift because we need domain [-Fs/2, Fs/2]
% ---------

% -------
% Now filter the signal
y = filter(H, 1, x);
Y = abs(fftshift(fft(y)));
% -------

That's it :)


In the Sampling Theorem, why are the image frequencies at n*fc not a problem?

The "alias" frequencies you see appearing above the Nyquist limit $f_s/2$ do not exist in a proper interpretation of the digitised signal; they are appearing due to incorrect interpolation of the samples during reconstruction. If you correctly interpolate the signal from the samples, those frequences will not be there.

Remember, a key condition of the Nyquist–Shannon sampling theorem is that it can perfectly record the input signal only so long as it contains no frequencies above $f_s/2$. A corollary of this is that the output signal you produce from correctly interpreting the samples must not contain such frequences either. The core of the theorem is that if you sample a signal under this condition, and then plot a sequence of points corresponding to each sample, there will be only one curve you can interpolate that both passes through every point and also contains no frequencies above $f_s/2$; that curve is the same as the original input.

The difficulties and details of adhering to this condition and properly interpolating the output are what cause problems like the one you're seeing.

In the case of your output, you're probably reconstructing the signal using piecewise interpolation, producing a "stair-step" output that looks like this:

(Image credit: https://wiki.xiph.org/Videos/Digital_Show_and_Tell)

That's certainly one way of doing it, but it (or at least it used alone) is the wrong way. Such an interpolation clearly has a lot of frequency information above $f_s/2$ (the alias frequencies) and is significantly different the original input signal meeting the conditions above.

It is, however, often an intermediate signal during the interpolation process; it turns out that one way of doing the interpolation fairly accurately is to to generate a signal using piecewise interpolation and then filter out all components of it above $f_s/2$, which will leave us with the original waveform.

You can do this in your first diagram by simply erasing everything above $f_s/2$, leaving only the original 5 Khz frequency you sampled.

For a more in-depth explanation of this, I recommend Monty Montgomery's, "D/A and A/D | Digital Show and Tell" video (also on YouTube). There's also a text version if you're short on time, but the video demonstration is signifantly more clear. I would consider this almost mandatory viewing for anybody involved in any details of digital sampling, audio or otherwise.


How to use a FIR on a signal with a drifting central frequency?

Without understanding more detailed results on what i(t) is, one solution is to estimate either the phase or frequency first, and then pass that result through a filter.

An approach to do the frequency estimation digitally directly on the I (real) and Q (imaginary) samples of x(t) is the cross product frequency discriminator. After hard limiting x(t) to remove the AM components (which would result in samples on the unit circle in the complex plane), the frequency can be estimated using:

$$F \propto  I_1Q_2-I_2Q_1$$

Where $I_1$,$Q_1$ and $I_2$,$Q_2$ represent successive samples. This is because the cross product is proportional to the phase between two vectors, so in this case the phase is the phase versus time, and therefore frequency. This is shown in the math and diagram below.

Frequency is the change in phase versus time ($d\phi/dt$). To determine the phase between two hard limited successive samples in time $V_1$ and $V_2$ given by:

$V_1 = e^{j\phi_1} = I_1+jQ_1$

$V_2 = e^{j\phi_2} = I_2+jQ_2$

For small angles, the phase $\phi_2-\phi_1$ can be reasonably approximated by the imaginary component of the conjugate product $(V^*_1) (V_2)$:

$V^*_1V_2 = e^{-j\phi_1}e^{j\phi_2} =e^{j(\phi_2-\phi_1)}$

For small angles $cos(\phi) \approx \phi$,

So therefore for small angles

$\phi_2-\phi_1 \approx Im(e^{-j\phi_1}e^{j\phi_2})$

$= Im((I_1-jQ_1)(I_2+jQ_2))$

$=I_1Q_2-I_2Q_1$

The above is extremely simple with low processing requirements useful for hardware implementation but depends on the small angle criteria for accurate angle estimation. This will work well as long as the sampling rate is high enough that the phase between samples doesn't deviate significantly. Alternatively if further processing is feasible, an ATAN2 operation can be instead be done on the two successive I,Q samples for full 4 quadrant accurate phase estimation (and therefore frequency given the phase difference between two successive samples over a much wider frequency range).

A simple delay and multiply frequency discriminator could also be used, followed by a  filtering of the resulting frequency versus time waveform. I have further detailed an example of using a simple delay and multiply frequency discriminator in this recent answer: Demodulation of FSK signal

Both approaches require hard limiting your waveform first to remove any amplitude variation of your signal. A differentiator as done through the differencing in these approaches is a high pass filter, the hard limiter removes all AM noise, and the post detection filtering will remove the noise induced "jumps" that you see in the result.


HINT:

$(-1)^n=e^{jn\pi}$. What is its DTFT?
Multiplication with $(-1)^n$ corresponds to convolution in the frequency domain with the DTFT you obtained in step 1.


How does sampling a frequency result in a duplicated spectrum?

if you consider the spectrum of analog signal as the figure in your picture.Sampling in time domain is same as multiplying a pulse train to  the analog signal.This means that the frequency response of pulse train is convoluted in frequency domain. The frequency response of pulse train is again a pulse train,hence the spectrum gets duplicated in frequency domain due to its convolution with the pulse train.


I think there are a couple things going on here.

First, the flatness of your baseline plot is a little deceiving because the tones that you have added (1 - 500 Hz) are all almost exact multiples of the fundamental frequency of the fft (fs/nfft = 1.024 Hz/bin).

If you increased nfft by a factor of 10, for example, the finer frequency resolution of the zero-padded FFT would show something similar to your New FFT Plot.  Here it is with nfft = 10240:

Because the tones are centered well on the FFT bins, there is no spectral leakage.  If you were to change the tones so that they were for example at multiples of 0.75 Hz, you'd see something like this:

It looks like a mess, but what's really happening is that each tones is getting smeared across several bins and the tones are constructively and destructively interfering.  Here's a single tone at 0.75 Hz with nfft=10240:

This spectral leakage occurs because the FFT assumes that its input is repeated forever, so if you have a sinusoid whose period is a multiple of the fundamental frequency of the FFT everything looks as you'd expect, but if not then there will be a discontinuity between the first sample and the last sample.  This effect can be mitigated using window functions: https://en.wikipedia.org/wiki/Window_function.

However, the main issue that you're seeing is the result of aliasing.  If you want to avoid aliasing, you need to low-pass filter your signal before downsampling.  Otherwise, any frequency content in the original signal that is above the Nyquist rate of your down-sampled signal will alias.  In your case, the sampling rate after down-sampling is 1024/10 = 102.4 Hz, so all of the tones you have above 51.2 Hz in your original signal will alias to frequencies below 51.2 Hz and interfere with the true signals at those frequencies.

Changing your down-sampling factor to 8 makes the effect easier to demonstrate because the aliased tones will land right on top of other true tones.  If your original signal only has tones below the 64 Hz, all tones are present in the down-sampled signal:

for f in np.arange(1, 64):
...
step = 8
x2 = x1[::step]   # down-sampling by a factor of 10??
fs2 = fs/step

Now let's add a tone at 65 Hz:

for f in np.arange(1, 66):

There is no longer any tone at 63 Hz or 65 Hz.  The 65 Hz tone has aliased and destructively interfered with the true tone at 63 Hz to cancel it out:

Now, let's change the phase of the original tone by 180 degrees:

for f in np.arange(1, 66):
x1 += np.sin(2*np.pi*f*t1 + np.pi/2)

There's still no tone at 65 Hz, but the one at 65 Hz is now constructively interfering, doubling the magnitude in the 63 Hz bin.

So if you extend this to 500 Hz, you can start to see what's going on.  The fix is to to apply a low-pass filter before you down-sample.  The combined steps of low-pass filtering and then down-sampling is called decimation.  Here's the scipy function to do it: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.decimate.html.


How does amplitude scale when FFT size is smaller then signal duration?

When you perform simple DFT of your signal with length $N$, then you have to scale it properly. Normally it is done by dividing magnitude by length of spectrum vector, which is $N$. Another thing is that normally we are interested in one-sided spectrum, therefore you need to multiply your above result by 2 to compensate energy. This will assure, that values of magnitude are representative.

In case when you are using windows different from rectangular (Hamming, Blackman, etc.) you are not dividing by $N$ - you must scale it by sum of your window samples.

Therefore when you are windowing your signal, you get a new signal with a length of a window. Thus in your case amplitude will be scaled down by 12.

Example for 1Hz sine wave and two window types is shown below. Please note, that frequency axis is chopped at 3Hz. You can see that amplitude of peaks is equal to unit amplitude of signal. Also in case of Hamming window side-lobes are much lower than in case of rectangular, sacrificing main-lobe width at the same time. This can be seen better in logarithmic scale.

MATLAB code used to produce result:

clc, clear all, close all
f0 = 1;        % fundamental frequency
fs = 50;       % sampling frequency
T = 6*pi;      % upper time limit
t = 0:1/fs:T;  % time vector
N = length(t); % number of samples

% 5 Hz sinusoid with unit amplitude
s = sin(2*pi*f0*t);
% Windowing signal
win = hamming(N)';
s_win = s.*win;
% Magnitude spectra
S = abs(fft(s));
S_win = abs(fft(s_win));
% Frequency vector
f = (0:fs/N:round(fs/2));
% Take one-sided spectrum
S_h = 2*S(1:length(f));
S_win_h = 2*S_win(1:length(f));
% Normalisation
S_h = S_h / length(s);
S_win_h = S_win_h / sum(win);
% Time domain signal plots
subplot(2,1,1)
plot(t,s)
xlim([0 T]); grid on; hold on;
plot(t,s_win,'r');
legend({'Rectangular window','Hamming window'})
xlabel('Time'); ylabel('Amplitude')
% Frequency domain plots
subplot(2,1,2)
plot(f, S_h)
grid on; hold on;
plot(f, S_win_h, 'r');
legend({'Rectangular window','Hamming window'})
xlabel('Frequency'); ylabel('Amplitude')
xlim([0 3])


What is the purpose of the filtering the image data in the PNG file format

The usage of a word filtering here is somewhat of a misnomer. The filtering operation is typically associated with waveform processing for various applications, as pulse shaping in the transmitter or noise suppression in the receiver. Image filters are also used, but primarily not for better compressibility of image files. In the context of your reference of PNG file format standard, it is possible that this operation would be more aptly named data preprocessing.
A Wikipedia article on PNG file format is more explicit on the compressibility of filtered PNG image data:

An image line filtered in this way is often more compressible than the
raw image line would be, especially if it is similar to the line
above, since the differences from prediction will generally be
clustered around 0, rather than spread over all possible image values.

To further clarify the effect of PNG data filtering on data compressibility, I cite the LZ77 algorithm article:

LZ77 algorithms achieve compression by replacing repeated occurrences
of data with references to a single copy of that data existing earlier
in the uncompressed data stream.

The less variability of "the differences from prediction" that are "clustered around 0", the more "repeated occurrences of data" exist "in the uncompressed data stream". Therefore, the more effective is the LZ77 algorithm's work.
The LZ77 is an algorithm the PNG standard uses for compression.
Notice also the statement from  W3C description of the PNG file format:

Filters are applied to bytes, not to pixels, regardless of the bit
depth or colour type of the image.

Again, the Wikipedia article is helpful to better understand this statement w.r.t. compressibility:

An image line filtered in this way is often more compressible than the
raw image line would be, especially if it is similar to the line
above, since the differences from prediction will generally be
clustered around 0, rather than spread over all possible image values.
This is particularly important in relating separate rows, since
DEFLATE has no understanding that an image is a 2D entity, and instead
just sees the image data as a stream of bytes.


Regarding your existing steps:

Do you have some overlap between adjacent blocks? It is common in feature extraction systems to have some overlap between adjacent blocks, so that a short transient event can be correctly captured (if it is right at the end of a block, it will be right in the middle of the next one).
Your idea of discarding blocks is dangerous. You can do this to roughly identify the beginning and end of a speech segment; but silence recognition can also be built into the recognition model.

Regarding the data size:

Your signal is split into frames. For each frame, you compute a 1-D MFCC vector. So in the end, you have a sequence of 1-D MFCC vectors. This is the data on which you train your HMMs.


This answer is about using sines instead of triangles in the bounds and bounce concept. There is no way to do it so that these defining conditions of the original concept are all satisfied:

While the boundaries stay constant the output is a pure sine wave (triangle for the original algorithm) that goes from the upper boundary to the lower boundary.
Moving the boundaries suddenly does not create a discontinuity in the output.
Moving the boundaries suddenly does not create a discontinuity in the derivative (slope) of the output.

Say the sine wave has reached the upper boundary and the boundaries are updated. Sine waves of other amplitudes cannot have zero derivative at that amplitude value, so the scheme fails.

But you can make it work if you constrain the boundary changes to times of zero crossings of the sine wave:

You just have to decide on the slope at zero crossing and your sine turns out as:  $$\text{boundary}\times\sin\left(\frac{\text{slope}\times\text{time}}{\text{boundary}}\right).$$


Using ROC to find stability of system in specific example

because even though it has more poles than zeros

WRONG. The Z-transform transfer function will always have equal number of poles and zeros. Your poles are at $z = -5/4$ and $z = 1/4$. Zeros are at $z = -1/2$ and $z = \infty$. Since ROC will always be concentric circle region without including poles, there are 3 possible ROC for the given transfer function.
$$
ROC_1 = 5/4 \lt |z| \lt \infty
$$
$$
ROC_2 = 1/4 \lt |z| \lt 5/4
$$
$$
ROC_3 = 0 \le |z| \lt 1/4
$$
$ROC_1$ represents an unstable but causal system (since it can include $z \rightarrow \infty$ so terms like $z^{-k} \,\,with\,\, k \ge 0$ will only be present (only negative or zero powers of $z$). And because only negative powers of $z$ are present, system is non-zero for $n \ge 0$. So it is causal. Example - $H(z) = a_0 + a_1z^{-1}+a_2z^{-2}+a_3z^{-3}+a_4z^{-4}+...$

$ROC_2$ represents stable but non-causal system because it includes unit circle $|z| = 1$ but does not include $z \rightarrow \infty$. Terms with $z$ raised to positive power can be there along with negative power of $z$. So you can rule out $z=0$ and $z \rightarrow \infty$.Example - $H(z) = a_0 + a_1z^{-1}+a_2z^{-2}+a_3z^{-3}+a_4z^{-4}+a_{-1}z^{1}+a_{-2}z^{2}+a_{-3}z^{3}$

$ROC_3$ represents an unstable and non-causal system because it does not include unit circle and terms with only positive power of $z$ are present (because ROC includes $z=0$ and does not include $z \rightarrow \infty$. Example - $H(z) = a_0 + a_{-1}z^{1}+a_{-2}z^{2}+a_{-3}z^{3}+...$

In your case, the ROC might be $ROC_3$.


Are complex exponentials the only eigenfunctions of LTI systems?

All eigenfunctions of an LTI system can be described in terms of complex exponentials, and complex exponentials form a complete basis of the signal space. However, if you have a system that is degenerate, meaning you have eigensubspaces of dimension >1, then the eigenvectors to the corresponding eigenvalue are all linear combination of vectors from the subspace. And linear combinations of complex exponentials of different frequencies are not complex exponentials anymore.

Very simple example: The identity operator 1 as an LTI system has the whole signal space as eigensubspace with eigenvalue 1. That implies ALL functions are eigenfunctions.


Ok, so the question is not "Can the Kalman gain be bigger than 1?" but "Can the Kalman gain for this system be bigger than 1?"

In that case, I suggest you look at the evolution of the Kalman gain given your system and see if there are any conditions where it does. Perhaps for this particular system it can't be. You'll need to specify all the parameters of your signal model to do that.

Of course it's possible!

Why do you think the KF gain is less than 1? What the Kalman gain is depends on the system model and the data being processed.

There are plenty of examples showing gains greater than 1. See, for example, the image below taken from here.


Laplace of step and integration are same?

This is because the impulse response of an integrator is $h(t)=u(t)$. The output which is the convolution with the impulse respoponse is
$$y(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau$$
and with $h(t)=u(t)$ it becomes
$$\begin{align}
y(t)&=\int_{-\infty}^{\infty}x(\tau)u(t-\tau)d\tau\\
&=\int_{-\infty}^{t}x(\tau)d\tau\tag{1}\end{align}$$
where $(1)$ is resulted from the fact that $$u(t-\tau)=\begin{cases}0& \forall \tau>t\\
1 & \text{otherwise}\end{cases}$$
The transfer function (which is the Laplace transform of impulse response) is $$H(s)=\mathcal{L}\{h(t)=u(t)\}=\frac{1}{s}$$


Why is FFT-based convolution efficient only for signals of large size?

The direct convolution approach has a complexity proportional to $n^2$
The FFT based approach has complexity $n \log(n)$.
Since there are unknown pre-factors and the complexity only holds asymptotically, the only thing one can infer from the complexities is that only for large enough $n$ the FFT approach becomes more efficient.
The $n=500$ threshold is thus probably determined empirically.


Typical parameters and stability of Chebyshev filters

I design Chebyshev first and second type filters

You are not going to see a lot of effect with this low order filters. If you want to study ripple and attenuation crank it up.

  What are typical values of maximum ripple in the passband and minimum attenuation in the stopband? I struggle to find optimal value of those parameters.

This really depends on your application. You can pick something easy for starters: 1 dB ripple and 40 dB attenuation. That's probably not useful for a real application but at least you can make some graphs

  I know that IIR filters are often unstable,

IIR filters are only unstable if the poles are on or outside the unit circle. Chebyshev are inherently stable unless you made a mistake. Implementation of a stable filter can indeed be unstable. That is heavily dependent on filter representation, topology, numerical accuracy and pole locations.  A 6th order filters in scipy should not be a problem unless you are trying to use fix point or the cutoff frequencies are extremely low.


This can be implemented, without needing to access the real and imaginary parts of $z$ individually, by writing

$$\operatorname{Re}(z) + \operatorname{Im}(z) = \frac{(1 + i)(z^* - iz)}{2},$$

where $i$ is the imaginary unit and $z^*$ is complex conjugation.

Factoring out a complex constant, $\xi = 1/2 + i/2$, turns this into a slightly tidier expression,

$$\operatorname{Re}(z) + \operatorname{Im}(z) = \xi(z^* - iz).$$

These expression can be derived by applying the trigonometric identities

$$\cos(x) = \operatorname{Re}(e^{ix}) = \frac{{e^{ix} + e^{-ix}}}{2}$$
$$\sin(x) = \operatorname{Im}(e^{ix}) = \frac{{e^{ix} - e^{-ix}}}{2i}.$$

Perhaps this will give you some inspiration for choosing a good function name.


Moving average frequency response over an image

OK.  There's a lot of things going on in your question, and they aren't really pertinent to the title.  Here goes:

... but that post uses a discrete time Fourier transform, whereas I've seen in a book the use of the discrete Fourier transform ...

OK.  The book was taking a short cut.  Hopefully it explained what it was doing and why.  For an infinitely-large image, the discrete-time Fourier transform (DTFT) is exactly correct; for any image the discrete Fourier transform (DFT) will have undesired edge effects.

the only difference is that the discrete Fourier transform version is obtained by sampling the discrete-time Fourier transform $H(\omega)$ in the points $\omega_k = \frac{2 \pi k}{N}$.

This is incorrect.  It is a not-bad approximation if some conditions are met:  For a time-domain signal $x(n)$ that is only nonzero for some interval $n \in [n_{min}, n_{max}]$ and $n < N$, and which has a DTFT whose spectrum is very small for $\omega \notin [-\pi, \pi]$, then the approximation is good to the extent that the spectrum outside of $[-\pi, \pi]$ is small.

I don't understand ... why the DFT has "meaning" only if $k \in [-\frac{N}{2}, \frac{N}{2}]$ and not $k \in [0, N-1)$.

The DFT has meaning for any contiguous interval of $N$ points (and even sets of $N$ points that aren't contiguous, but meet certain conditions -- but that's getting really weird).  But that meaning gets further and further away from easy physical truths the further you get from $k \in [-\frac{N}{2}, \frac{N}{2}]$
For any DFT, $H(k) = H(k \mod N)$.  This drops out of the defining math.  It is why the DFT only has "meaning" for some set of $N$ points -- that's just all there is.  This is a reflection of the fact that the DFT is a lossless transform, and you start with $N$ points.  Once you look at $N$ output points, there cannot be more information in additional points.
For any DFT of a real-valued (not complex) signal, $H(k) = H(-k)^*$, where $X^*$ is the complex conjugate of $X$.  So for that case, while the DFT output has "meaning" outside the interval $k \in [-\frac{N}{2}, \frac{N}{2}]$, it's not a terribly useful meaning.  In fact, because $H$ is complex, and it contains $N$ non-zero real or imaginary components in the interval $k \in [0, \frac{N}{2}]$, that's often the only part you look at if you're starting with real-valued data (and there's FFT routines that are designed to take real-valued data and be just a bit faster than their fully-complex counterparts because of the inherent symmetry).

moreover, how does the size of the moving average filter relate to the frequency response?

Generally speaking, if you have two filters of the same overall shape, the broader one has a narrower frequency response, and visa-versa.  There's an exact relationship when you're talking about a continuous-time signal and taking it's Fourier transforms.  Everything becomes hand-wavy when you get into discrete time, but the principle more or less applies.
So -- a wider moving average filter has a lower cutoff frequency, pretty much proportional to the number of points.


let's first explain what warp is: if you apply LK for two images and you get say u=2 and v=3 for a certain pixel, in this case applying warping of one image is to increase the x-coordinate of that pixel by 2 and increase it's y-coordinate by 3, and then make this for all other pixels in the image using the associated u and v.

In step two you need to warp one image using the u and v motion vectors that you get from the first step, if your estimate to u and v from the first step is accurate then after applying the second step the two images will coincides on each other, and this means that if you apply Lucas Kanade algorithm on these two images(step 3) you will get u=0 and v=0 because the two images coincides on each others. but in real situations you will not get accurate u and v from the first iteration. so after applying the warping there will be a smaller motion between the two images and if you apply Lucas Kanade algorithm between these two images you will get a smaller u and v and you refine the motion vectors by the following:

$u_(next) = u_(previuos) + u_(current)$

$v_(next) = v_(previuos) + v_(current)$

the algorithm terminates when $u_(current)$ and $v_(current)$ are zeros or very small.

this is iterative LK without pyramids.

EDIT :

I get the idea from the following code which implements iterative LK with pyramids explained in the following paper:

B.D. Lucas and T. Kanade, "An Iterative Image Registration technique, with an Application to Stero Vision," Int'l Joint Conference Artifical Intelligence, pp. 121-130, 1981.


How does cascading low pass modules affect pass band ripples?

The passband region of an equi-ripple filter is the region where its magnitude is in the interval $[1-\delta,1+\delta]$, where $\delta$ is the maximum approximation error. At the passband edge the value of the magnitude must equal $1-\delta$ (as shown in your drawing). The passband edge is the last frequency where the magnitude is still inside the interval $[1-\delta,1+\delta]$, and since the magnitude then falls off towards the stopband, the magnitude at the passband edge must be at the lower limit of the passband magnitude, which equals $1-\delta$.

The total magnitude response of a cascade of $M$ identical filters is

$$|H_{total}(\omega)|=|H(\omega)|^M$$

Consequently, the total magnitude at the passband edge is $(1-\delta)^M$. This value is then equal to $1-\delta_{total}$, which gives

$$\delta_{total}=1-(1-\delta)^M$$


This can be useful to you: Oscar Dalmau Cedeño, Mariano Rivera, Teresa Alarcón: Bayesian Scheme for Interactive Colourization, Recolourization and Image/Video Editing. Comput. Graph. Forum 29(8): 2372-2386, 2010.  ISSN: 0167-7055


reference parseval theorem limited signal

I don't think you will find a reference for this. Your are simply stating that "integrating over a whole bunch of zeros results in 0". In other words, if
$$x(t) = 0, x\not\in [0,T]$$
then
$$ \int_{-\infty}^{+\infty}x(t) dt = \int_{-\infty}^{0}x(t) dt + \int_{0}^{T}x(t) dt + \int_{T}^{+\infty}x(t)dt = \int_{0}^{T}x(t) dt $$
since
$$ \int_{-\infty}^{0}x(t) dt = \int_{-\infty}^{0}0 dt=0$$ and

$$\int_{T}^{+\infty}x(t)dt = \int_{T}^{+\infty}0dt = 0$$

If that doesn't work, than I don't know what would


Frequency modulation audio synthesis: clicks/pops when changing modulator frequency dynamically

If you change modulator frequency, do note that it changes a sudden jump in phase, as it changes the result of the equation where the phase is calculated from frequency. This is largely unnecessary, as having a certain frequency means the phase change is constant per sample. Thus, if you change the frequency to be a phase increment per sample, then changing the frequency does not cause a sudden phase jump but just changes the phase increment from that point on. That's how FM synth chips work.


Cepstral windowing is the most obvious way to extract a spectral envelope, you have a description here.

  I've looked at cepstral enveloping but its very fiddly and I can't
  seem to get it to give the right results.

In the linked description, try lowering the lowpass "cut-off" sample (nc) until you retain the amount of detail you want. You can get useful implementation examples here.

  I've also tried LPC but again it eems to miss a lot of the minor peaks
  while enveloping the major peaks perfectly.

Isn't that the point of extracting the spectral envelope? By LPC are you referring to the Linear Prediction Spectral Envelope?


Transfer function tries to describe a system using the Fourier Transform (More generally the Laplace Transform / Z Transform).
Yet in order for the Transform function to dependent on the frequency and amplitude only the input system must be LTI system.
Hence there is no Transfer Function in the Fourier Domain to your system which can be described the way you want.
Linear Time Invariant (LTI) System
By Convolution Theorem a convolution of LTI System can be described by a transfer function in Fourier Domain.
What do mean by a transfer function? We mean having $ H \left( \omega \right) $ which is independent of the input.
As for any output and input we can calculate the Fourier Transform and have the ratio of them as the Transfer Function. But only for LTI the result of this ratio will be constant.


Tracking object that is missing in some frames

The short answer is: use your optical flow algorithm to get the best-fit warping of some number of previous images into your current image.  Then assume that your object is moving on a line through the images, and do some kind of robust best-fit to that line.  When the object reappears it is more likely to be near the predicted projection of the line than not.

You didn't give enough specifics in your question to really say more than that.

I assume by "LK" you mean "the Lucas-Kanade optical flow algorithm?"  But what are you really doing?  Are you taking image pyramids of each image and then doing multiple iterations of warping at each level of the pyramid?  Are you really doing the simplest form of Lucas-Kanade or are you actually doing one of the varieties of Kanade-Lucas-Tomasi where you are concentrating on "corners"?

What is your description of an "object"?  Is it really just a "feature" (a corner point) or a small set of features?

You may need to use a more accurate and robust optical flow algorithm.  For example, some variation that does dense estimation (rather than the sparse estimation done by Lucas-Kanade) and that uses more robust distance metrics than least squares.

For example:

Secrets of optical flow estimation and their principles
Sun, D., Roth, S., and Black, M. J.,
IEEE Conf. on Computer Vision and Pattern Recog., CVPR, June 2010.


Meaning of concentric circles on the spectrogram

In Understanding FFT Overlap Processing, p. 10, it is suggested that uniform of constant windows, for certain hop patterns, with repeated pulses, could generate such artifacts.


Maybe this thread helps you in this Q - Prewarping both resonant frequency $f_0$ and bandwidth (or $Q$) when using bilinear transform
H(s) is the s-plane transfer function

Analog prototype of RBJ peak filter (Octave):

K  = 10^(gain/40);
w0 = 2*pi*fc;
% H(s) = (s^2 + (K w0/Q)*s + w0^2) / (s^2 + w0 s/(K*Q) + w0^2)
BC = [1 K*w0/Q w0^2];
AC = [1 w0/(K*Q) w0^2];


How to get a real signal back to a non-complex form after processing it in frequency domain (using fft)

First I would advise against your filtering approach, which is the "Frequency Sampling" method of filtering (if you are using multiplying your FFT by a target response directly) which has poor performance and efficiency. (The frequency sampling approach will provide your exact solution at your FFT bin centers only, and then a lot more ripple in between vs the algorithms I suggest below.)  However if you need to be in the frequency domain anyway, and are multiplying by the FFT of your desired filter coefficients, this would be ok as long as you are properly dealing with the circular convolution involved.

Also I would use a real filter unless you are intentionally trying to get an asymmetric spectrum response (meaning positive and negative spectrums are differrent, which implies a complex signal-- given that is what you are trying to avoid then there is no reason to use a complex filter).

To implement a real filter, with a real signal and a real output, consider using either the Parks-McClellan or Least Squares algorithms, with design tools readily available in Matlab, Octave and Python, and then implement your filter as an FIR filter in the time domain.

If you need to implement filtering in the frequency domain, you can force it to be a real filter by ensuring that the filter is conjugate symmetric; for an FFT the 0 bin is the center (DC value), and then the samples from 1 to N/2-1 should be conjugate symmetric to the samples from N/2 to N-1; where sample 1 is the conjugate of sample N-1, sample 2 is conjugate of sample N-2 etc...


Replicating Research: Feature Vectors of Audio Samples

Here's my interpretation, and I'll summarize using Python's Librosa module:

Each frame is n_fft=512 samples (23.2 ms at sr=22050 Hz sampling rate), and hop_length=256 samples to get 50% overlap. Let y be the audio signal vector. When you run

mfcc_mat = librosa.feature.mfcc(y=y, sr=sr, n_fft=512, hop_length=256, n_mfcc=40)

you'll get a matrix with 40 rows ("coefficients") and $\approx$len(y)/hop_length columns ("time slices"). The authors only use the first 25 rows, and collapse the time dimension in 11 different ways (min, max, median, mean, etc.). Here's the (incomplete) code snippet that generates the 275-length feature vector:

metric_min = np.min( mfcc_mat[0:25,:], axis=1)
metric_max = np.max( mfcc_mat[0:25,:], axis=1)
metric_median = np.median( mfcc_mat[0:25,:], axis=1)
...
feature_vector = np.hstack((metric_min, metric_max, metric_median,...))

etc.

https://librosa.github.io/librosa/generated/librosa.feature.mfcc.html
https://librosa.github.io/librosa/generated/librosa.feature.melspectrogram.html#librosa.feature.melspectrogram


If H is representation of shift invariant operator then its a convolution operator. If you are implementing your algorithm in MATLAB, there is a intrinsic function with name convmtx2 which exactly do what you want.

Algorithms which warp image don't use a linear operator (a matrix) to do the job (considering large size of this operator it's not practical, even it was practical it's not an efficient method). They usually find the coordinates of warped image which generally don't match the pixel coordinates of original image then obtain pixel values using some interpolation method like nearest neighbor or bi-linear to obtain the pixel value from its neighbors.

If you want to do it by using a linear operator you have to choose an interpolation method then find the transformed coordinates and compare them to the original coordinates, then find the pixel weights regarding the interpolation method and store these weights in their corresponding location in the matrix (in nearest neighbor method give weight 1 for the nearest pixel and zero for others, in bi-linear method 4 non-zero weights for 4 pixel which are around the transformed coordinate and zero for other pixels).

Good luck.


Beginner - use FFT to display spectrum/equalizer

I think we need a little more information to fully answer your question.  Most FFT routing (like FFTW) provide many different ways to call the transform, and the return result will depend on exactly what routine you call, and how you call it.

Having said the standard disclaimer, if you perform an FFT on $N$ real data points, you will get back $N/2$ unique complex points.  The basic algorithm though will assume your $N$ real points are $N$ complex points with zero imaginary parts.  Many routines will thus return $N$ complex points, but the second half will be (in the case of real input) hermitian symmetric (that is, point $N/2+i$ will be the complex conjugate of point $N/2-j$).

So, start by recognizing that your result is complex numbers, not real, and you will need to at least plot the "square" ($xx^*$) or each point.  If you have the real and imaginary points returned as $X_r$ and $X_i$, then you want to plot $X_r^2+X_i^2$.

So you are probably pretty close.  You will want to consider a few other things: windowing your input data, and averaging your output data.  To do the first, you need to perform the transform on weighted data, rather than your direct data.  See http://en.wikipedia.org/wiki/Window_function for more info.  You will simply hand the FFT function data of the form $W_k   {x_k}$ instead of just the $x_k$ (where $k$ is the index number).  This will be particularly important if you use a fractional frequency in your sin function (try it, you'll see).

The averaging step involves averaging the squared output. You will simply want to average a few outputs, point by point, together.  You will also often want to plot the output as a logarithm though, instead of simply plotting the squared output.


If I understand correctly, you are trying to create a chirp an arbitrary instantaneous frequency, $f(t)$.
In that case, the phase of the chirp is given by $\varphi(t) = 2 \pi \int f(t)dt$ and your signal in time domain is given by $s(t) = \sin(\varphi(t))$.
For example, in case of a linear chirp  you mentioned, $f(t) = f_0 + \beta t$ (where $\beta$ is the sweep rate) and therefore $s(t) = \sin \left(2 \pi (f_0 t+\frac{1}{2}\beta t^2) \right)$.


What went wrong?

How do you know it's wrong? Could be that you have way more energy in the first harmonic than at the fundamental.
Debug rule #1: test the code first with a signal where you already know the answer (say a mixture of a few sine waves). Once you have made sure that your code produces the expected result EXACTLY you can move on to real world signals.
Your frequency resolution seems insufficient. You only have 10 Hz of resolution, that will not give a lot of detail at 30 Hz. Keep in mind that the energy at each particular frequency will be "smeared out" across a few neighboring FFT bins, so you want enough resolution to see that.
Your lowpass seems way to steep. Together with filtfilt you get a 20th order filter which will have enormous amount of time domain ringing which filtfilt() unfortunately truncates and the truncation can create significant errors


A PID that can handle 2 inputs for my application

If you require two outputs, which belong to a movement in x and the other in y direction, you could start with the assumption, that both movements are independent an simply implement two PID controller, each caring for one direction.

It gets interesting at the point where you translate the output of these two controller to signals to the actuators. Here it may be wrong to assume that both movements are independent. If you give more insight on that, we can give better ideas for control algorithms

(more a comment then an answer, but I cant comment)


FPGA is only user to verify the initial RTL which is specifically written to do each of the task listed below. Once the verification is done there is a tape out and everything is implemented in hardware accelerators as below. FPGA is wasteful of area and too power hungry which cannot be afforded in modem (low power device)

The LTE channel estimation is done as follows:

The FFT of the incoming samples at receiver is taken post cyclic
prefix removal. This is done in an FFT kernel in Hardware (specifically written RTL, not in  FPGA)
Then this result is fed to a memory where all the    subcarriers are
stored for this particular symbol (other symbols are    stored in
this memory as well).
Then firstly all positions at which    CRS    (Cell Reference Symbol)
is available are extracted for the symbol, the channel is estimated (MMSE or others) at the subcarriers corresponding to the CRS. This is done in a hardware accelerator. See figure below

Then there is an interpolation done in between the CRS subcarriers to
find the the channel estimate at non CRS subcarriers (wiener
filters), this is called interpolation in frequency. This is also
done in Hardware accelerators (i.e. RTL is written specifically for
this)
Now, CRS (per antenna port) is only transmitted in 2 symbols in a
slot(for most of the transmission modes). It  is transmitted in the
first and $5^{th}$ symbol of a slot for normal CP.
All the symbols for a full subframe  are buffered and frequency interpolation is done at all symbols in which CRS is transmitted as explained above (see figure below)

Then there is an interpolation done in time to get channel estimates
for all symbols in between the symbols containing CRS, to get a
channel estimate for all the subcarriers in the entire time frequency grid. This is also done in the Hardware accelerators
The CRS references positions are at every 6th subcarrier in frequency
(linked to the coherence bandwidth) and similarly in time linked to
coherence time

Once there are channel estimates available for all subcarriers  of 1 full subframe. It is sent for equalization/demodulation (this is also done in Hardware accelerators)

The reason for using Hardware accelerators over doing in Firmware is speed. The reason for doing in accelerators over FPGA is area and power. accelerators have a smaller footprint and more power efficient

To the first question, by itself there is no limitation that an LTE channel can have for the phase shift. But there are limitations as to what can be estimated unambiguously by the receiver. Consider this situation: The maximum frequency offset that an LTE system can estimate is +/-7.5KHz (because subcarrier spacing is 15 Khz) anything above this frequency offset will wrap around to within the above specified range. Now since frequency offset shows up as phase in the signal multiplied with channel estimate the ambiguity will remain in unwrapped phase.

As to how then we measure frequency offsets more than +/-7.5Khz, that's another story :)


1-) The value at t = 0 of an ideal step function converges to $\frac{f(0-) + f(0+)}{2}$ at 0 if you approach the question from a Fourier transform perspective (see Fourier transform of a discontinuous function).

2-) The analogous function in discrete time is still a unit step defined as:
$$
u[n] =
\begin{cases}
0 & \text{if n <  0} \\
1 & \text{if n >= 0} \\
\end{cases}$$

The concept of infinite frequency does not convey to discrete time representation.  In discrete time, frequency is represented by sequences of complex exponentials rather than continuous complex exponential functions.  Complex exponential sequences are not unique for frequencies above 1/2 your discrete time sampling rate so your frequency range is limited to $[0, \frac{f_s}{2}]$, where $f_s$ is your sampling rate.

See sampling theory, the nyquist frequency and a discussion of the discrete Fourier transform for a more detailed discussion.


How to decide the number of cepstral coefficients?

From CMUSphinx FAQ:

  There are various types of MFCC which differ by number of parameters,
  but not really different for accuracy (it might be a few percent worse
  or better).

  The interpretation of MFCC (Roughtly introduced Alan V. Oppenheim and
  Ronald W. Schafer. From Frequency to Quefrency: A History of the
  Cepstrum. IEEE SIGNAL PROCESSING MAGAZINE) is not applicable as such,
  and the use of 12 or 13 coefficients seem to be due to historical
  reasons in many of the reported cases. The choice of the number of
  MFCCs to include in an ASR system is largely empirical. To understand
  why any specific number of cepstral coefficients is used, you could do
  worse than look at very early (pre-HMM) papers. When using DTW using
  Euclidean or even Mahalanobis distances, it quickly became apparent
  that the very high cepstral coefficients were not helpful for
  recognition, and to a lesser extent, neither were the very low ones.
  The most common solution was to “lifter” the MFCCs - i.e. apply a
  weighting function to them to emphasise the mid-range coefficients.
  These liftering functions were “optimised” by a number of researchers,
  but they almost always ended up being close to zero by the time you
  got to the 12th coefficient.

  In practice, the optimal number of coefficients depends on the
  quantity of training data, the details of the training algorithm (in
  particular how well the PDFs can be modelled as the dimensionality of
  the feature space increases), the number of Gaussian mixtures in the
  HMMs, the speaker and background noise characteristics, and sometimes
  the available computing resources.


SNR calculation in audio reconstruction and what is an acceptable value

The common practice in the source separation community is to report 3 metrics defined in Vincent et al's paper, called SDR (signal to distortion ratio), SIR (signal to interferences ratio), and SAR (signal to artifacts ratio).

The reason for this is that in source separation applications, there can be two types of "noise" in the reconstructed signal: noise due to mis-separation, called interferences (for example, if you want to extract the vocals from a song, this would be a residual of the background music), and noise due to the reconstruction algorithm itself, called artifacts (for example, glitches due to the STFT phase estimation process). The former is sometimes referred to as "crosstalk" in other contexts.

There is usually an application-dependent trade-off between the two - for example in a Karaoke application (vocal removal), it is not a big problem if the audio accompaniment has still a few traces of the original voice (low SIR), as long as the result is pleasing to the ears and glitch-free (high SAR). Most separation algorithms have "knobs" that will allow you to balance SAR and SIR. For example, increasing the STFT window size tends to increase SIR but decreases SAR.

The estimation of these metrics is described in the paper by Vincent et al, and implemented in the bss_eval toolbox. Of course, access to the original individual sources you want to recover is necessary! You will have to evaluate your algorithm on a data-set for which the "right answer" (original individual track) is known.

Regarding the acceptable value of the metrics, you need to compare your results with those published in the academia for the same problem you are attempting to solve. Keep in mind that these are objective metrics, and as such, do not take into account the perceptual quality of the noise. It might be possible that an algorithm with a -12dB SAR will be more pleasing to the ears than an algorithm with a -24dB SAR if the former produces artifacts well spread in time and below the masking threshold of the reconstructed signal. SDR/SIR/SAR can be used in development to rule out bad ideas and tune your algorithms, but for a commercial application, listening tests are ultimately necessary.


You can create a sine wave in numpy.  Think of it as if you were sampling data from a signal, only you are sampling a mathematical function:

sample_rate = 1024
dt = 1.0/sample_rate
t = np.arange(sample_rate)*dt  # 1 second of samples
freq = 5
amp = 1.0
sine1 = amp*np.sin(2*np.pi*freq*t)
sine2 = .5*np.sin(2*np.pi*15*t)
sinsum = sine1 + sine2


Find a stable transfer function $G(z)$ such that $|G(z)| = |H(z)|$

Your idea of replacing the pole $p$ outside the unit circle (i.e., $|p|>1$) with a pole $1/p^*$ inside the unit circle is correct. However, as mentioned in this answer, you need to be careful with the scaling. In order to obtain the correct scaling, a factor $(z-p)$ in the denominator must be replaced by a factor $\pm |p|(z-1/p^*)$. In that way the magnitude of the transfer function will be unchanged on the unit circle, i.e., for $|z|=1$. This should also be added to the problem formulation:

  ... such that $|G(z)|=|H(z)|$ for $|z|=1$.

Clearly, the new stable transfer function $G(z)$ can be written as the original transfer function $H(z)$ cascaded with an allpass filter $A(z)$, since $|A(z)|=1$ for $|z|=1$:

$$G(z)=H(z)A(z)\quad\Longrightarrow\quad |G(z)|=|H(z)|,\quad |z|=1\tag{1}$$

For the given problem, the transfer function $A(z)$ is given by

$$A(z)=\frac{z-3}{1-3z}\tag{2}$$


How do you handle negative frequencies in the power spectrum of a complex signal?

For a real signal, content at the negative frequencies generated by using the DFT is redundant. This is due to the well-known property of real signals with respect to the Fourier transform family: their transforms are Hermitian symmetric. That is, for any real signal $x[n]$,

$$
X[k] = \sum_{n=0}^{N-1}x[n] e^{-j2\pi n k / N} = (\sum_{n=0}^{N-1}x[n] e^{-j(-2\pi n k / N)})^* = X[-k] = X[N - k]
$$

So, if your input is a real signal, all of the information is in the positive frequency bins; the negative frequencies can be discarded for many applications. However, there is no such proprety for general complex signals. They can have power spectra that are asymmetric about zero frequency, so you may not discard any of the frequency bins of a complex signal's power spectrum without loss of information.


What are the audio equivalents of images like "Lena", "Mandrill", and "Cameraman"?

The closest example I can think of is the beginning of Suzanne Vega's "Tom's Diner" which has been used for the mpeg-1 layer 3 development, and is still occasionally used to demo audio codecs.


Does Wigner-Ville distribution provide the highest resolution?

The statement is true if you restrict it to all TF distributions in the Cohen Class. That is rather easy to see, because the Wigner distribution generates all other distribution in that class by convolution with a non-negative time frequency kernel and therefore attains the minimal uncertainty.

Beyond the Cohen class things get much more complicated. First of all, the very simple concept of uncertainty in the Cohen class is not necessarily useful for non-bilinear distributions. So it's not obvious how to even generalise the question. So statements that attempt to be more general than the Cohen class first need some clarification about mathematical terminology and possibly even concepts. Outside the Cohen class you will also find distributions without time-frequency translation invariance. If they're bilinear they can have a TF uncertainty that is just as good as the Wigner-Ville distribution. Dropping bilinearity too, you can even beat WV with reassignment distributions of all sorts.

So summing up, the statement is true with the (possibly implcit) restriction to the Cohen class, which is the most common theoretical framework for TF distributions. Without that restriction it is not accurate.


Produce musical signal from instantaneous frequency

The f[i] * i together determine the phase and frequency at every i so sine wave is only rendered correctly if there is no speed change, so the algorithm must be changed to allow for continuous time i and any change in frequency must not cause a phase jump.
So try a phase accumulator approach, where phase for each i is incremented as required by the frequency.


If X plus Y is longer than twice your sample rate, 20 seconds, but not an exact rational multiple of your sample rate, then with a large enough number of samples you might be able to make a frequency and duty cycle estimate by calculating the statistics and spectrum of the sample vector.


How can we prove the correctness of the integration property of the Laplace transform?

You are right that the argument in the proof is not correct, or at least misleading. The fact that $e^{-st}$ becomes zero as $t\to\infty$ is true for any $s$ with $\textrm{Re}\{s\}>0$. However, the integral
$$\int_0^{\infty}f(\tau)d\tau$$
might not exist, so we can't just compute the limit by claiming that the first term becomes zero.
We have to compute
$$\lim_{t\to\infty}-\frac{e^{-st}}{s}\int_0^{t}f(\tau)d\tau\tag{1}$$
The limit $(1)$ only equals zero for $s$ inside the region of convergence (ROC) of $f(t)$, i.e., we don't just require that $e^{-st}\to 0$ for $t\to\infty$, but we require it to decay sufficiently fast to compensate for the growth of the integral.


Why autocorrelation can be more efficiently calculated using the fft

The cross-correlation between two functions $f(t)$ and $g(t)$ can be seen as a convolution of $f(t)$ and $g(-t)$. The auto-correlation is of course a special case where $f=g$. This operation in the frequency domain corresponds to a multiplication. Therefore, an efficient way to compute the cross-correlation is to compute the Fourier Transform of both $g(t)$ and $f(t)$, multiply them in the frequency domain and the anti-transform the result, to go back to the time domain.

This approach apparently looks more computationally expensive, however in the discrete case, if we use the fast algorithm of the FFT (Fast Fourier Transform) to compute the DFT (Discrete Fourier Transform), which uses $N\log N$ arithmetic operations, it leads to more efficient algorithm than the actual computation of the cross-correlation function in the time domain, which requires $N^2$ operations, where $N$ is the number of samples. So the number of operations is much lower if $N$ is particularly high, which is often the case ...

You can watch this video if you want to know more about this https://www.youtube.com/watch?v=iTMn0Kt18tg


PSD estimate using FFT - scaling and units

Have a look at this paper, section 9 goes over how to scale or normalize the FFT output. For the PSD you need to divide by the effective noise bandwidth (ENBW) which is

$$
ENBW = f_s\frac{S_2}{S_1^{2}}\\
S_1 = \sum_{j=0}^{N-1}w_j\\
S_2 = \sum_{j=0}^{N-1}w_j^{2}\\
where~w_j~is~the~jth~sample~of~the~window~function
$$


An energy or zero-crossing criterion might help in detecting segments containing no voice, but it won't help in removing the noise superimposed with speech. Since you are mentioning removing noise for improving speech analysis, this doesn't look like the right technique to use.

If you compute your energy using small windows, down to summing a few adjacent coefficient, you will see fluctuation in energy which are just due to the periodic variation of the waveform itself! As a rule of thumb, the window size on which you compute an energy feature must be at least one magnitude bigger than the period of the signal to analyze. Assuming a f0 of 200 Hz, this gives you a minimal integration window of 50ms, or 400 samples at 8kHz.


Image registration / fusion optimal similarity metrics

I'm not sure what is meant here by the "optimal" similarity metric for a given situation. The textbook you linked to references the following thesis:

  [87] P. A. Viola, Alignment by maximization of mutual information. Ph.D. thesis, Massachusetts Institute of Technology, 1995.

...but I can't find it anywhere online after a quick search.

From my point of view:

Sum of square differences measures the difference in pixel intensities, so it is useful when pixel intensities are very similar in both images
Normalized cross-correlation measures linear correlation between pixel intensities, so it is useful when image intensities are linearly correlated
Mutual information is a measure of general statistical dependence between pixel intensities, so it is useful when image intensities have a non-linear statistical dependence

If SSD works for a given pair of images, I don't see why NCC or MI wouldn't do the job just as well, since they are more general. They may be slower since the metric takes longer to compute, but I don't understand why one would be more "optimal" than the other.


Why FFT of a rectangular pulse has imaginary part while the sinc function is real?

The even part of a real signal transforms to the real part of the DFT.  The odd part of a real signal transforms to the imaginary part of the DFT.  For example, the DFT of cosine, an even function, is strictly real, while the DFT of a sine, an odd function, is strictly imaginary.


Purely theoretical question about idelal filters and infinite oscillations

Note: historically the Gibbs phenomenon was described over a continuous-time domain periodic square wave which was constructed from a finite number of its Fourier series coefficients, however, one can interchange the domains and observe the same behaviour on the frequency response of a finite length discrete-time filter with its impulse response samples designating the Fourier series coefficients.

This question is essentially asking the consequences of the unavoidable Gibbs phenomenon on the mathematical construction of an ideal brick-wall filter whose frequency domain definition apparently (and confusingly) does not involve the expected Gibbs ripples...?

The practical answer is that as the filter length approaches infinity the peak magnitude of the ripples does not reduce, but its duration shrinks to a single point at the edges of the passband. From Riemannian integration point of view, this will not have an effect on the filter's output as if there's no Gibbs phenomenon.

On any finite length approximation of an ideal brickwall filter, however, the ripples do exist, and will alter the actual output from that of an ideal filter's computation. The better solution to reduce those ripples is to make use of a window instead of increasing the filter impulse response truncation length.


How to implement Signal transmitted by RF source

I assume this would be as follows:

s*exp(j*2*pi*f*n)

Where:

n: sample index
s: baseband signal, magnitude and phase for each sample
f: carrier frequency in cycles/sample

Note how the carrier as a discrete signal is given in cycles per sample. This is the normalized frequency and is related to the desired frequency $F$ in Hz by dividing by the sampling rate $f_s$ as follows:

$f = F/f_s$


Unfortunately I don't have Matlab now but I'm thinking your problem has to do with linspace. Remember linspace gives the actual first and last values always. When you're doing the FFT binning I don't believe you want the boundaries to be pinned that way. I'd try using the colon operator instead and appropriate starting and delta values.

EDIT:
If you look at the definition of the DFT it is:
$$X(k)=\sum_{n=0}^{N−1}x(n)e^{j\frac{2πkn}{N}}$$

This shows that the first DFT output at $k=0$ is at baseband. Every subsequent frequency is at an interval of $F_s/N$. So if you had 10 points the bin centers would be from $0...(\frac{N-1}{N}F_s=0.9F_s)$. For 1000 points it would span $0...0.999F_s$. Linspace, on the other hand always leaves the first and last points exact. You need the start and spacing exact, or to adjust the ending point with linspace.


Converting poles/zeros to differential/difference equation solutions

Okay, I'm gonna attempt an answer to this while I am still convinced the question is ill-formed.  I am answering the title of the question without really referring to the text of the question, which I don't understand.

First, for discrete-time:  A general IIR (or even FIR in which all of the poles are at the origin, where $z=0$) transfer function is:

$$\begin{align}
 H(z) \triangleq \frac{Y(z)}{X(z)} & = H_0 \frac{(z-q_1)(z-q_2)(z-q_3)...(z-q_M)}{(z-p_1)(z-p_2)(z-p_3)...(z-p_N)} \\
\\
 & = H_0 \frac{z^M (1-q_1 z^{-1})(1-q_2 z^{-1})(1-q_3 z^{-1})...(1-q_M z^{-1})}{z^N (1-p_1 z^{-1})(1-p_2 z^{-1})(1-p_3 z^{-1})...(1-p_N z^{-1})}  \\
\\
 & = z^{M-N} \ \frac{H_0 \ (1-q_1 z^{-1})(1-q_2 z^{-1})(1-q_3 z^{-1})...(1-q_M z^{-1})}{(1-p_1 z^{-1})(1-p_2 z^{-1})(1-p_3 z^{-1})...(1-p_N z^{-1})}  \\
\\
 & = z^{M-N} \ \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + ... + b_M z^{-M} }{ 1 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + ... +  + a_N z^{-N} }
\end{align}$$

$N$ is the number of poles, $p_i$, and $M$ is the number of zeros, $q_i$.

The difference equation comes from moving a few of the above things around.

$$\begin{align} Y(z) & \left(  1 + a_1 z^{-1} + a_2 z^{-2} + ... + a_N z^{-N}  \right) \\
& \quad = X(z) z^{M-N} \left( b_0 + b_1 z^{-1} + b_2 z^{-2} + ... + b_M z^{-M} \right) \end{align}$$

$$\begin{align} Y(z) = \ & b_0 z^{M-N} X(z) + b_1 z^{M-N-1} X(z) + b_2 z^{M-N-2} X(z) + ...  \\
 & \quad ... + b_M z^{-N} X(z) - a_1 z^{-1}Y(z) - a_2 z^{-2}Y(z) - ... - a_N z^{-N}Y(z) \end{align}$$

Inverting the Z Transform, this becomes

$$\begin{align} y[n] = \ & b_0 \ x[n+M-N] + b_1 \ x[n+M-N-1] + b_2 \ x[n+M-N-2] + ...  \\
 & \quad ... + b_M \ x[n-N] - a_1 \ y[n-1] - a_2 \ y[n-2] - ... - a_N \ y[n-N] \end{align}$$

For a causal system or filter, we cannot look into the future for any input sample $x[n]$.  That means $M \le N$ so that there is no negative delays on the $x[n+M-N]$ term.

Often $M=N$ (the numbers of poles and zeros are the same) and you get the more familiar

$$\begin{align} y[n] = \ & b_0 \ x[n] + b_1 \ x[n-1] + b_2 \ x[n-2] + ... + b_N \ x[n-N] \\
 & \quad - a_1 \ y[n-1] - a_2 \ y[n-2] - ... - a_N \ y[n-N] \end{align}$$

and the more familiar transfer function

$$\begin{align}
 H(z) & = H_0 \frac{(z-q_1)(z-q_2)(z-q_3)...(z-q_N)}{(z-p_1)(z-p_2)(z-p_3)...(z-p_N)} \\
\\
 & = \frac{H_0 \ (1-q_1 z^{-1})(1-q_2 z^{-1})(1-q_3 z^{-1})...(1-q_N z^{-1})}{(1-p_1 z^{-1})(1-p_2 z^{-1})(1-p_3 z^{-1})...(1-p_N z^{-1})}  \\
\\
 & = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + ... + b_N z^{-N} }{ 1 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + ... +  + a_N z^{-N} }
\end{align}$$

That demonstrates where the difference equation (and the feedforward, $b_i$, and feedback, $-a_i$, coefficients) come from the poles, $p_i$, and zeros, $q_i$.

The problem remaining, which is the hard problem, is that I didn't bother to multiply it out.  So you have to blast this out:

$$\begin{align}
 H_0 \ & (1-q_1 z^{-1})(1-q_2 z^{-1})(1-q_3 z^{-1})...(1-q_M z^{-1})  \\
 & \quad = b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + ... + b_M z^{-M}
\end{align}$$

and

$$\begin{align}
  & (1-p_1 z^{-1})(1-p_2 z^{-1})(1-p_3 z^{-1})...(1-p_N z^{-1})  \\
 & \quad = 1 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + ... + a_N z^{-N}
\end{align}$$

You can see right away that $b_0 = H_0$ and that $a_0=1$ (I never really had an "$a_0$" term anyway), but you have to blast out the multiplication of the binomials to get the polynomials.  You can also see that $b_M = (-1)^M H_0 q_1 q_2 q_3 ... q_M$ and $a_N = (-1)^N p_1 p_2 p_3 ... p_N$   That's not a simple expression, but I hope you understand what you must do here.


You may not need to delve into scipy or opencv.  The following solution uses pil and numpy:

from PIL import Image
import numpy

#define a greyscale threshold which we will use later
threshold = 230

#open image as greyscale with PIL and convert to numpy array
with Image.open('your_image.png').convert('L') as img:
    arr_gry = numpy.array(img)

#create 1-D array of max column values
colmax = numpy.zeros(arr_gry.shape[0])
colmax = numpy.min(arr, axis=0)

#analyse the bands array and define start and end of each band
band_limits = {}
band_list = []
prev_column = 0
for i,column in enumerate(colmax):
    if column > threshold and ('start' not in band_limits):
        band_limits['start'] = i
    if column <=threshold and prev_column > threshold:
        band_limits['end'] = i
        band_list.append(band_limits)
        band_limits = {}
    prev_column = column

###The following is for illustration purposes only

#convert to RGB array
arr_clr = numpy.stack([arr_gry,arr_gry,arr_gry], axis=2)

#Colour the bands and mid points
for band in band_list:
    arr_clr[:,band['start']:band['end']] = [255,0,0]
    arr_clr[:,(band['start']+band['end'])//2] = [0,0,255]

#Convert array to Image
out = Image.fromarray(arr_clr)
out.show()

This produces the following image

Obviously you can use opencv if desired to open, view and save the images.  But in this example all the processing is done in numpy so pil is more than capable.


What are the measurable factors of image sharpness?

Some references on image sharpness metrics:

Encoding Visual Sensitivity by MaxPol Convolution Filters for Image Sharpness Assessment, IEEE Transactions on Image Processing, 2019
A Fast Approach for No-Reference Image Sharpness Assessment Based on Maximum Local Variation, IEEE Signal Processing Letters, 2014
Image Sharpness Assessment Based on Local Phase Coherence, IEEE Transactions on Image Processing, 2013
$S_3$: A Spectral and Spatial Measure of Local Perceived Sharpness in Natural Images,  IEEE Transactions on Image Processing, 2012


It's exactly the same method, only the second is converted to (power) dB ($\text{dB}:\,x\mapsto 10\log_{10}x$):

Let's convert your first formula to dB:

$$\begin{align}
P_{1,\text{dB}} &= 10\log_{10}P_1\\
&=10\log_{10}\left(\frac{\Delta t}{N}\lvert\text{FFT}\rvert^2\right)\\
&=10\log_{10}\left(\frac{\lvert\text{FFT}\rvert^2}{f_sN}\right)\\
&=10\log_{10}\left(\frac{\lvert\text{FFT}\rvert^2}{\sqrt{f_sN}^2}\right)\\
&=10\log_{10}\left(\left(\frac{\lvert\text{FFT}\rvert}{\sqrt{f_sN}}\right)^2\right)\\
&=10\cdot2\cdot\log_{10}\left(\frac{\lvert\text{FFT}\rvert}{\sqrt{f_sN}}\right)\\
&=20\log_{10}\left(\frac{\lvert\text{FFT}\rvert}{\sqrt{f_sN}}\right)\\
\end{align}$$

The $+30$ just happen because your seconds source is considering $\text{dBm}$ (decibel of milliwatt) instead of $\text{dBW}$ (decibel of watt), assuming your input values are in Volt over a impedance of 1Ω or so.

Note: you nowhere in your question stated the physical units of what your input data is; be very careful when "implying" something can be represented as $\text{dBm}$ or $\text{dBW}$. That makes a statement about your signals origin (namely, that your ADC or wherever the signal came from is calibrated, or will be converted to analog with calibrated equipment) if you're dealing with real-world observations/signals. This assumption is usually wrong – if you just got a good microphone and a good sound card (or a good receiver and a good SDR device) your digital amplitudes will be proportional, but not identical, to the voltages observed, or even more significantly, to the sound intensity (or EM field strength) that caused the electrical signal in the first place.


edge detection evaluation medical test images MATLAB

The closest resource I have come across is the BSD500 database:

http://www.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/resources.html

It is a set of natural images that have had the contours drawn by humans. When you download the test images there are also MATLAB files with all the contours.

However, the question is always "what is an edge?" Is it a step in brightness in the local image patch? Is a line an edge? Is it a boundary between semantic objects in a 2D image? The latter is what the above set tests. The best performing approaches for it also consider an abrupt change in texture to be an edge.

I guess the point is the performance of your edge detector is dependent on the application. One could argue that every edge detector is perfect as it perfectly detects edges according to its own definition. But that is just being cheeky.

You could also test some other properties, such as

Is it invariant to rotation? Does the same line at 0 degrees and X degrees give the same result?
Is it invariant to scaling? I.e. zooming of the image
Is it invariant to view point? I.e. warping of the image
Is it illumination invariant? Does it detect true weak edges as well as strong edges?

I also have a new edge detector that I've been working on. One of recent ideas I have had is to use the Mikolajczyk test set. It consists of eight sets of images with various transformations, including the above. Typically the test set is used to evaluate interest point descriptors.

The process would be

Detect edges in each image in the set.
Using the homographies, transform each detection image back to the first image in the set so they can be compared
Compare the edges between images using one of the measures you mention.

Ideally, if the edge detector is invariant to the various translations, the same edges should be found. This approach might have already been done in the literature but I haven't looked.


What are some promising areas of research in digital signal processing?

Choosing a novel topic is a bold and wise move. Your question requires a bit of a bet, yet it is interesting to me, as I have not thought about this for a while. I will first drop two dual thoughts (they could be updated):

shed light on black-boxes: the present successes of machine learning and deep networks shall not make us forger that signals are our (DSP community) sensing door to a globally unknown world. Truly understanding what we can understand with data and models remains at our core, and we can be much better at that (to make it more interpretable, robust, model-based, explainable, reliable etc.)
save electrons in devices and computations: sensing and computing are becoming pervasive (IoT, etc.), and this requires materials (from the extraction of rare-earth elements for sensors and computing devices) and resources (form energy-harvesting sensing to bit-wise or qubit-related fast algorithms). Saving power (efficient filters, compression, fast algorithms, etc.) are in our DNA.

DSP is a blended field of knowledge, where people of mixed origins (computing, acquisition, information theory, physics, optimization, computer vision, statistics) meet. I hope that we DSP people can innovate here.
Additional stuff, as they come:

The Future of Signal Processing: conference program (2017)


Since you indicated that the power spectrum of your background noise is flat, I'll assume it is white. A major drawback with your current approach is that you are discarding a large amount of the signal power; even with the effect of front-end bandlimiting shown by in your diagram by the exponential-rise step response, a single ADC sample near the end of the rounded pulse provides a snapshot of the receiver input that is rather localized in time. You can take advantage of more of the signal power by sampling at a higher rate and applying a matched filter at the higher sample rate.

Theory:

You can look at this as a relatively simple problem in detection theory. In each symbol interval, your receiver needs to decide between two hypotheses:

$$
\begin{eqnarray*}
H_0 &:& signal \ is \ not \ present \\
H_1 &:& signal \ is \ present
\end{eqnarray*}
$$

This sort of problem is often solved using Bayesian decision rules, which attempt to make the optimum decision according to some specified measure of risk. This provides a framework where one can optimally make detection decisions based on a flexible set of criteria. For example, if there is a large penalty to your system for failing to detect the signal if it is in fact present (i.e. you pick $H_0$ when $H_1$ is true), then you can build that into your decision rule if needed.

For a detection problem such as yours, where you're trying to decide between zeros and ones at the receiver output, the penalty is typically assumed to be equal (outputting a zero when a one was transmitted, and vice versa, "hurt equally"). The Bayesian approach in that case reduces to a maximum-likelihood estimator (also described here): you pick the hypothesis that is most likely, given the observation that your receiver makes. That is, if the quantity that your receiver observes is $x$, then it would generate a decision based on the hypothesis that has the largest likelihood function value. For the binary decision case, the likelihood ratio can be used instead:

$$
\Lambda(x) = \frac{P(x\ |\ H_0 \ is \ true)}{P(x\ |\ H_1 \ is \ true)} = \frac{P(x\ |\ signal \ is \ not \ present)}{P(x\ |\ signal \ is \ present)}
$$

Using the above model, for each observation of the channel $x$, the optimal receiver would decide that the signal was not present (therefore outputting a zero) if the likelihood ratio $\Lambda(x)$ is greater than one (and therefore the signal was most likely to be not present based on the observation), and vice versa.

What remains is a model for the signal of interest and any other components in the receiver detection statistic $x$ that could affect its decisions. For a digital communications like this, it might be modeled as follows:

$$
\begin{eqnarray*}
H_0 &:& x = N \\
H_1 &:& x = s + N
\end{eqnarray*}
$$

where $n$ is a random variable taken from some distribution (often assumed to be zero-mean Gaussian) and $s$ is a deterministic component of the observation that is due to the signal that you're looking for. The distribution of the receiver observable $x$, therefore, varies depending upon whether hypothesis $H_0$ or $H_1$ is true. In order to evaluate the likelihood ratio, you need a model for what those distributions are. For the Gaussian case referenced above, the mathematics look like this:

$$
\Lambda(x) = \frac{P(x\ |\ H_0 \ is \ true)}{P(x\ |\ H_1 \ is \ true)} = \frac{P(x\ |\ x = N)}{P(x\ |\ x = s + N)}
$$

$$
\Lambda(x) = \frac{P(x\ |\ H_0 \ is \ true)}{P(x\ |\ H_1 \ is \ true)} = \frac{e^{\frac{-x^2}{2 \sigma^2}}}{e^{\frac{-(x - s)^2}{2 \sigma^2}}}
$$

where $\sigma^2$ is the variance of the Gaussian noise term. Note that the additive signal component only has the function of shifting the mean of the resulting Gaussian distribution of $x$. The log-likelihood ratio can be used to get rid of the exponentials:

$$
\ln(\Lambda(x)) = \ln\left(\frac{e^{\frac{-x^2}{2 \sigma^2}}}{e^{\frac{-(x - s)^2}{2 \sigma^2}}}\right) = \left(\frac{-x^2}{2 \sigma^2}\right) - \left(\frac{-(x - s)^2}{2 \sigma^2}\right)
$$

Recall that our decision rule picked $H_0$ if the likelihood ratio was greater than one. The equivalent log-likelihood decision rule is to pick $H_0$ if the log-likelihood is greater than zero. Some algebra shows that the decision rule reduces to:

$$
\begin{eqnarray*}
x < \frac{s}{2} \rightarrow choose\ H_0 \\
x > \frac{s}{2} \rightarrow choose\ H_1
\end{eqnarray*}
$$

Note that if $x = \frac{s}{2}$, then both hypotheses are equally likely, and you would need to just pick one; this isn't a practical concern for continuously-valued signals, however. So, given a known signal amplitude $s$, we can detect its presence against a background of Gaussian noise optimally by setting a threshold $T = \frac{s}{2}$; if the observed value $x$ is greater than $T$, we declare the signal present and emit a one, and vice versa.

Practice:

There are a few practical issues that creep into this simple, toy theoretical example. One: just mapping the scenario that you described into a deceptively simple-looking model might not seem straightforward. Secondly, it's very rare that you would know the amplitude $s$ of the signal that you're looking for, so threshold selection requires some thought.

As I referenced before, noise is often assumed to be Gaussian because the normal distribution is so easy to work with: the sum of a bunch of independent Gaussians is still Gaussian, and their mean and variances just add also. Also, the first- and second-order statistics of the distribution are enough to completely characterize them (given the mean and variance of a Gaussian distribution, you can write its pdf). So, hopefully that's a decent approximation at least for your application.

There are two ways to improve the performance of the detector given the model described above: you can increase $s$ (i.e. increase the signal power), making it stand out more against the noise. You could decrease $N$ (i.e. reduce the amount of noise), reducing the amount of interference that makes the presence of $s$ unclear. Or, equivalently, you can think of the signal to noise ratio instead. To see its importance, let's go back to the theory for a second. What is the probability of a bit error given our decision rule?

$$
\begin{eqnarray*}
P_e &=& P(choose \ H_0 \ |\ H_1\ true) P(H_1\ true) + P(choose \ H_1 \ |\ H_0\ true) P(H_0\ true) \\
&=& \frac{1}{2} P(x < \frac{s}{2} \ |\ x = s + N) + \frac{1}{2} P(x > \frac{s}{2} \ |\ x = N) \\
&=& \frac{1}{2} F_{x\ |\ x = s + N}\left(\frac{s}{2}\right) + \frac{1}{2} \left(1 - F_{x\ |\ x = N}\left(\frac{s}{2}\right)\right)
\end{eqnarray*}
$$

where $F_{x\ |\ x = s + N}(z)$ is the cumulative distribution function of the distribution of the observation $x$, given that $x = s + N$ (and likewise for the other function). Substituting in the cdf for the Gaussian distribution, we get:

$$
\begin{eqnarray*}
P_e &=& \frac{1}{2} \left(1 - Q\left(\frac{\frac{s}{2} - s}{\sigma}\right)\right) + \frac{1}{2}  Q\left(\frac{\frac{s}{2}}{\sigma}\right) \\
&=& \frac{1}{2}  + \frac{1}{2} \left(-Q\left(\frac{\frac{s}{2} - s}{\sigma}\right) + Q\left(\frac{\frac{s}{2}}{\sigma}\right)\right) \\
&=& \frac{1}{2} + \frac{1}{2} \left(-Q\left(\frac{-s}{2\sigma}\right) + Q\left(\frac{s}{2\sigma}\right)\right) \\
&=& \frac{1}{2} + \frac{1}{2} \left(-Q\left(\frac{-SNR}{2}\right) + Q\left(\frac{SNR}{2}\right)\right) \\
&=& Q\left(\frac{SNR}{2}\right)
\end{eqnarray*}
$$

where $Q(x)$ is the Q function:

$$
Q(x) = \frac{1}{\sqrt{2 \pi}} \int_{x}^{\infty} e^{\frac{-z^2}{2}} dz
$$

(i.e. the tail integral of the standard normal distribution's pdf, or $1$ minus the distribution's cdf) and $SNR$ is the signal-to-noise ratio $\frac{s}{\sigma}$. The above function is a strictly decreasing function of $SNR$; as you increase the ratio of the signal amplitude $s$ to the noise standard deviation $\sigma$, the probability of making a bit decision error decreases. So, it behooves you to do whatever you can to increase this ratio.

Remember our assumption that the noise was white and Gaussian? That can help us now. If the noise is white and Gaussian, then the noise components contained in each observation are jointly independent of one another. An important property of independent random variables is that when you sum them together, their means and variances sum. So, let's consider another simple case, where instead of taking one sample per symbol interval, you take two, then sum them together. I'll assume for simplicity that the pulse shape is rectangular (not an exponential rise), so the signal component $s$ in each observation $x_1$ and $x_2$ is the same. What is the difference in signal to noise ratio between just a single observation $x_1$ and the sum of two independent ones?

$$
SNR_1 = \frac{s}{\sigma}
$$

$$
SNR_2 = \frac{2s}{\sqrt{2 \sigma}} = \sqrt{2}SNR_1
$$

So, the signal to noise ratio in the combined observation is larger than using only a single sample (under the assumption of equal signal component and equal-variance white Gaussian noise in both samples that we took). This is a basic observation that points out the potential benefits of taking more than one sample per symbol interval and integrating them together (which, for a rectangular pulse, is a matched filter). In general, you want to cover the entire symbol interval with samples so that your receiver "ingests" as much of the transmitted energy for each symbol, thus maximizing the SNR in the combined output. The ratio of symbol energy to the background noise variance $\frac{E_s}{N_0}$ is often used as a figure of merit when evaluating digital communications system performance.

More rigorously, it can be shown that a matched filter has an impulse response that is identical in shape (that is, "matched", with the only subtle exception being that the impulse response is reversed in time) to the pulse shape that the receiver sees (so it weights more strongly samples that have larger signal components). That shape is a function of the transmitted pulse shape as well as any effects induced by the channel or receiver front end, such as bandlimiting or multipath.

To implement this sort of arrangement in practice, you would convolve the stream of samples taken by your ADC with the time-reversed expected pulse shape. This has the effect of calculating the cross-correlation between the pulse shape and the received signal for all possible time offsets. Your implementation is aided by the precise time synchronization that you have available, so you'll know exactly which matched filter output samples correspond to correct sampling instants. The filter outputs at those times are used as the detection statistic $x$ in the theoretical model above.

I referred to threshold selection before, which can be a complicated topic, and there are many different ways that you can choose one, depending upon your system's structure. Selecting a threshold for an on-off-keyed system is complicated by the likely-unknown signal amplitude $s$; other signal constellations, like antipodal signaling (e.g. binary phase shift keying, or BPSK) have a more obvious threshold choice (for BPSK, the best threshold is zero for equally-likely data).

One simple implementation of a threshold selector for OOK might calculate the mean of many observations. Assuming that zeros and ones are equally likely, the expected value of the resulting random variable is half of the signal amplitude, which is the threshold that you seek. Performing this operation over a sliding window can allow you to be somewhat adaptive to varying background conditions.

Note that this is only intended to be a high-level introduction to the issues inherent in digital communications with respect to detection theory. It can be a very complicated topic, with a lot of statistics involved; I tried to make it somewhat easy to understand while keeping true to the underlying theory. For a better explanation, go get a good textbook, like Sklar's.


Detect quasi-parallel lines with spurs in a binary image

I would start with a pruning algorithm (see skeleton pruning) in order to erase small segments, then use:

the Hough transform, but segments are not really straights.
if the pruning was a success, then you separate each segment and use the first and last point to compute the orientation, which will define when line are parallel. You will have to introduce a tolerance because you want quasi-parallel lines.


Because $s-s_{\alpha1} = 0 $ and $s-s_{\alpha2} = 0 $ so poles: $s_1,s_2 = s_{\alpha1},s_{\alpha2} $
What numerical values you have for $s_{\alpha1},s_{\alpha2} $ are the What numerical values you are looking for.


Should we divide the power at 0Hz by two for Single-Sided Power Spectrum with FFT?

The factor of "2" comes from the contribution of the negative frequencies, which doesn't apply to DC.
A better way to plot this is to NOT multiply with 2 but plot the entire FFT range from -400Hz to +400Hz. There you will see three components: 1 W at DC and 0.25W each at -50Hz and +50Hz which is exactly what's happening here.
If you want to plot power spectrum one sided, you need to multiply by two EXCEPT for DC and Nyquist, which have no negative frequency equivalent.


Down sample incoming stream of data to plot a dynamic graph

Downsampling

Let's define our downsampling factor :
$$r = \frac{F_s'}{F_s}$$

Where $F_s$ is the original sampling frequency and $F_s'$ the reduced one.

If $F_s$ is a multiple of $F_s'$ (i.e. $1/r$ is an integer), then the dowsampling is straightforward : you just take one sample out of $1/r$. For example, if you want to half your sampling frequency, just discard every $2^{nd}$ sample and that's all (well not really, see next paragraph).

However, if that's not the case, you're going to have to interpolate between the samples. To make it simple, interpolating means trying to guess the value of a sample by knowing the value of the surrounding samples. You brain does interpolation all the time : when you see several points forming a line you instinctively draw the line in your head. Several types of interpolation exist, the more simple beeing the less precise. It usually depends on your application.

But there's more...

You've maybe heard of Shannon-Nyquist sampling theorem, which states that to sample a signal, you have to have a sampling frequency at least twice as big as your signal's highest harmonic. If you don't, it's going to cause something called aliasing. I won't explain this theorem here because there is plenty of information available everywhere.

But keeping this theorem in mind, you can probably understand that downsampling your signal might cause aliasing if you're not careful. The simple solution to that is to filter the signal prior to downsampling, at a cutoff frequency of $\frac{F_s'}{2}$. The steeper the filter's transition band, the less aliasing.

Back to your application

From what I understood, you have a signal that you want to map to a screen. For that, you need to downsample the signal to your screen's horizontal resolution. Let's call this resolution $H_x$ and your signal's sampling frequency $F_s$. For the sake of the example, let's fix $H_x=240\text{px}$ and $F_s=5\text{kHz}$.

You also have a time parameter (let's call it $t_r$) which sets the "zooming" of your scope, ie how long of a signal corresponds the length of the screen. Let's set it to $t_r=100\text{ms}$ (Disclaimer : I have absolutely no idea about the frequency of ECG so this is completely arbitrary).

That gives you :
$$
F_s' = \frac{H_x}{t_r} = \frac{240}{0.100}=2400\text{Hz}
$$
Note that that is not a multiple of $F_s$.

Applying what I said before, you have :
$$ r = \frac{F_s'}{F_s} = \frac{2400}{5000} = 0.48 $$

So first, you filter your signal with a brickwall at $f_c=\frac{F_s'}{2} = 1.2\text{kHz}$, then you downsample it by a factor of $0.48$ using interpolation.

Another thing to consider is the vertical resolution. For that it is a bit different as you are not downsampling the signal but quantizing it. However, the process is similar and I think you can extrapolate without our help.


How to distinguish between the different frequency domains?

This answer is primarily opinion. It seems to me that the most sensible thing to do is to make the transform of interest clear in your text. Then you can use "frequency domain" exclusively and unambiguously.

For Laplace transforms, "s domain" will work as well, and I do believe they are interchangeable. I do not believe "Fourier domain" would be correct usage in this case. "s domain" may be preferred if you wish to emphasize that real values of $s$ have a clear aperiodic interpretation (exponentials in the time domain).

If the transform is not yet clear, you could say "Fourier domain", I suppose. In that case the assumed spectrum would be the complex (possibly conjugate symmetric) form.

Still, if it is within your control, be explicit about the chosen transform.