Effectively extracting “real” data from a noisy dataset

Judging from your description alone, I'd say an algorithm to find appropriate 3-clusters of timestamps would be something like:

compare the three first elements of your sorted lists. Pick the earliest of the three.
Remove the element from its list.
determine difference of that element to the other two first

If both difference are below a "simultaneity" threshold, register a new "common event".

return to 1

Often, you'd also remove the two "matches" from the other two list in 3., but that reduces your ability to detect ongoing physical excitation of your sensor as one by post-processing your "common events" list.
Seeing you're only dealing with about a billion timestamps in total, this simplified variation of your distance-limited cluster detector should allow for quick experimentation¹.

But:
Reading through your question, it seems likely that one physical phenomenon can lead to multiple adjacent timestamps in a single list, so that would be another factor that might help to characterize real events vs false alarm. In that case, better modelling might be highly desirable, but as a rough idea:
Imagine your three lists as three parallel lines in three-dimensional space. They are arranged with equal distance to each other, i.e. there's a fourth parallel line in their "middle" that has the same distance to each one of them.
Now, imagine your timestamps are points on the respective lines.
For every point on your magical middle line, you draw a sphere around that point and count the time stamps in that sphere. Setting an appropriate sphere radius and threshold allows for good detection.
If you have a mathematical model for how observed timestamps relate to real phenomena, that can change the shape of your sphere, and make it "fuzzy", i.e. weighing points further away less than nearby.
In the end, you'd have a function that says "for each time on my middle line, this is the likelihood of a real event having happened"; you can then build a maximum likelihood detector.

¹ hint: This might be a bad stereotype, but you sound like a matlab user. If there's one easy thing matlab is slow at, it's doing actual m-code loops and removing heads from lists; it's always preferable to just work at the end of a contiguous vector of data, since all you need to do to "remove" an element is change your own notion of length.


time frequency localization using wavelet transform

What information are you trying to extract from your signal?

  I want to know the frequency (heart rate) at every time point

If this is the information you want then any sort of frequency analysis is unlikely to be very useful. It will show you that you have a 1 or 2 hz periodic signal with a particular frequency profile, but you already know all of this information so that won't be especially illuminating.

You are probably better off with peak detection, or local minima/maxima detection to extrapolate beats.

If you can provide a sample of your data, and/or more details about the exact information you are trying to extract then it will help us to find a more precise solution to your problem.


How do I warp a 1 dimensional signal in frequency domain?

The mapping in this specific figure is basically a "lookup table" that maps $\omega$ to $\tilde{\omega}$ and the process in detail is outlined in section 5 of the paper that you link.

But on the point of applying that to a series, a function that effects such a transformation that is controlled by a parameter $\alpha$ is $f(x) = x^\alpha, x\in \left[0..1 \right]$ and $x \in \mathbb{R}$.

For example (Using Octave but easily extendable to other platforms):

function y=f(x,a) y=x.^a;end;
x = 0:0.001:1;
plot(x,f(x,1),x,f(x,2),x,f(x,1./2));
xlabel('x');ylabel('f(x)');
legend("x","x^2","x^{0.5}");
grid on;

Would give you:

Notice here that this is slightly different than the transform that you have with respect to the behaviour of the function with $\alpha$. I am using this for illustrative purposes here.

Notice also that because $x \in \left[0..1 \right]$, we can transpose that function to whatever limits we like. In your particular example it is $\pi$ or the number of bins of your Frequency Transform. Let's call that $N_{FT}$.

To apply the transform, you can either discretise the function or evaluate it on the fly. For example:

x_lookup = round(N_FT .* f(x/N_FT,alpha));

In this case, your input frequency transform bin would be x and that would have to now appear in bin x_lookup.

Notice here that x/N_FT to make sure that x stays within the range $0..1$ (assuming of course its lowest value is 0 here).

Again, section 5 contains all the details for specific f(x) functions you will need to effect this kind of mapping for a specific application.

Hope this helps.


Discrete-time, second-order section, direct-form II filter

This is Direct Form 1 but it's still cascaded SOS.  it's in C, but you should be able to translate it to Python easily.  the register keyword is not important.
typedef struct
{
float filterCoefficients[5];    /* b0, b1, b2, a1, a2   (a0 is normalized to 1) */
} sectionData;

void filterBuffer(int32 Nsamples, int Nsections, float* y, float* x, float* filterStates, sectionData** sections)
{
sectionData** section_ptr;
float sampleValue, state1, state2, *state_ptr, *coef_ptr;

for (int n=0; n<Nsamples; n++)
    {
    section_ptr = sections;                     // reset section pointer
    state_ptr = filterStates;                   // reset state pointer

    sampleValue = x[n];                         // get input sample

    state2 = *state_ptr++;                      //   x[n-2]
    state1 = *state_ptr--;                      //   x[n-1]
    *state_ptr++ = state1;                      //   x[n-1] -> x[n-2]
    *state_ptr++ = sampleValue;                 //   x[n]   -> x[n-1]

    int i = Nsections;
    while (--i >= 0)
        {
        coef_ptr = &((*section_ptr++)->filterCoefficients[0]);  // point to section filter coefficients

        sampleValue =  *coef_ptr++ * sampleValue;   //   b0*x[n]
        sampleValue += *coef_ptr++ * state1;        //   b1*x[n-1]
        sampleValue += *coef_ptr++ * state2;        //   b2*x[n-2]

        state2 = *state_ptr++;                      //   y[n-2]
        state1 = *state_ptr--;                      //   y[n-1]
        *state_ptr++ = state1;                      //   y[n-1] -> y[n-2]

        sampleValue += *coef_ptr++ * state1;        //  -a1*y[n-1]
        sampleValue += *coef_ptr++ * state2;        //  -a2*y[n-2]

        *state_ptr++ = sampleValue;                 //   y[n]   -> y[n-1]

        //     at this point sampleValue is the output of the section just completed
        //     and will be the input to the section about to be computed unless this
        //     is the last section, in which case sampleValue is the output.
        //     state1 and state2 are the output states of the section just completed
        //     and will be the input states of of the section about to be computed.
        }

    y[n] = sampleValue;                             // store output sample
    }
}


Thanks for the answer. You are correct - the training sequence is known by the time you establish any communication.

In case if somebody stumbles upon this in the future, here's how to get this done:

The synchronization burst contains an extended training sequence. It's longer than the regular one and there is only one sequence. You are supposed to use that to equalize, demodulate and decode the SCH burst. Among other useful information, it contains the color code of the BTS, which identifies the training sequence used by that BTS.

The procedure is nicely explained in this book:

Sajal Kumar Das: "Mobile Handset Design"


Missing delay in heavyside step function

We know that $\displaystyle
\frac{\sin(\theta (n+1))}{\sin(\theta)} u[n+1]$ has value $0$ for $n < -1$ since $u[n+1]=0$ for $n < -1$. At $n=-1$, $u[n+1]$ jumps to value $1$, but $\sin(\theta (n+1))\bigr|_{n=-1} = \sin(\theta (-1+1))$ has value $0$, and so $\displaystyle
\frac{\sin(\theta (n+1))}{\sin(\theta)} u[n+1]$ just happens to equal $\displaystyle
\frac{\sin(\theta (n+1))}{\sin(\theta)} u[n]$.

No need to convolve with $\delta[n+1]$ (or subtract $\delta[n+1]$ as Ch. Siva Ram Kishore suggests, why I don't know).


Speech recognition using mfcc and neural network

Usually your training data and your test data are different, otherwise your results will be unrealistically good. However, for a first test of your implementation, you can simply use your training data as test material. In this case, the recognition performance should be very good. If this is not the case you know that something went wrong.


use wavelet transform to extract frequencies from given signal

A simple (hanning-windowed) FFT reveals four peaks that stand out far beyond the noise. This appears to be what you are looking for.


Problem with applying elliptic filter to ecg signal

Numerical problems. 64-bit double precision floating point is not nearly enough to implement  this filter in the way you have done it.
Your filter is extremely steep: The poles are way too close to the unit circle and the order is way too high to implement it in transfer function form (which is a bad anyway).
Things to do

Implement the filter as cascaded second order section
Review the requirements for your filter. It seems way steeper than it needs to be. Such an aggressive filter will destroy a lot of time domain details since the filter has massive amounts of time domain ringing. The settling time alone is multiple seconds and by using filtfilt you double it again & create an enormous non-causality.

Filter design requires a lot of trade offs, make sure you understand them and optimize to the specific requirements of your application.


Add echo as described by difference equation to audio signal (in MATLAB)

The impulse response $h[n]$ can be obtained by inspection:
$$
h[n] = \delta[n] + k \delta[n - d]
$$

In MATLAB notation this would be:

h = [1, zeros(1,d-1), k];

The output is computed by:

y = conv(x, h);

Note that the time sample indices for the vector y will be as follows, assuming that both x and h are indexed where the first element corresponds to time $n=0$:

ny = 0:(length(x) + length(h) - 2);


What's the impact of aliasing in the time domain?

I read Footnote 8,

  "The complication is that because of the sampling, the total system
  is not time-translation invariant and so does not have a unique
  ‘impulse response’ – the response is slightly different according to
  the position of an original impulse relative to the sampling points."

as a cumbersome way of saying that aliasing corrputs an LTI system such that its response can no longer be trusted to be linear and time-invariant.

If I were to rephrase things informally, I'd say that:

When a signal is undersampled in the frequency domain or has inadequate frequency resolution (aka freq domain aliasing), its spectrum has overlapping tails.

That means that it is no longer possible to correctly reconstruct the time-domain signal. The spectrum tails do not go to zero, but are folded back. This inversion of the tails are what we know as aliasing and manifest themselves in the time-domain signal as impulse responses (a short-duration time-domain signal) that mix with the original signal.

These mixtures are perceived as distortion/roughness/etc. Their constructive/destructive effect give you your "quasi-random" amplitude modulations.

The magnitude & position of aliasing in the time domain is related to the degree of frequency domain undersampling. See https://ccrma.stanford.edu/~jos/sasp/Example_2_Time_Domain.html

Finally, regarding "time smearing", the impulsive (quickly rising) parts of a waveform contain many high-frequency harmonics that when undersampled in the frequency domain, fold back as described above. However, I wouldn't call this distortion smearing, in the transient-smearing, or pre-ringing sense.

Hope that helps!


2D Fourier Transform of Rotated Discrete Domain Signal

The Rotation Property means that the choice of coordinate direction would not affect the spectrum of the signal itself.

To perform the rotation, I think it is better to do in the spatial domain since frequency domain is a complex one and it will just increase the complexity.


How do I estimate the derivative of the current position, when I have only values from past to present?

Is the derivative between x[n] and x[n-1] useful enough, hence x'[n] = x[n] - x[n-1] ? Or do you need a better estimate?


Why do we discard imaginary part of the phase spectrum?

Phase spectrum by itself is a real function of frequency and, therefore, does not contain any imaginary parts. It is the fourier transform which is in general a complex function of frequency, henceforth, may include real and imaginary parts. But neither the magnitude spectrum nor the phase spectrum contain(s) any imaginary parts.

Phase spectrum is a real function. Whatever methods/functions (correct or not) you use to obtain it does not and shall not change this fact. $H(e^{j\omega}) = |H(e^{j\omega})| e^{j\phi(w)} $ where the $\phi(w)$ contains the phase spectrum and the absolute value part is the magnitude spectrum. Both are real. You shall use a proper function and method in your environment to obtain those signals from $H(e^{j\omega})$


Should I use ideal or non-ideal filters for offline filtering?

For most practical applications, filter responses that approach the behavior of the ideal "brick-wall" filter are overkill. I know it's tempting to try to design a really sharp filter when you've got all the time in the world (i.e. for offline applications), but if you really look at the characteristics of your problem, you can most likely get away with something much more reasonable.

Another good reason not to use a brick-wall filter: its impulse response is infinite in length, and it has the form of a $sinc$ function. When you apply such a filter to your signal, you may notice long-duration ringing in the signal at the output of the filter. This comes from the fact that the $sinc$ function in the filter's impulse response is infinitely long and doesn't decay very fast; the resulting effect is most likely not desirable. In general, filters with long frequency responses are not very suitable for analysis of short signals, as the filter output is dominated by the transient behavior of the filter.

The aforementioned effects aren't limited strictly to ideal filter responses; if you design a filter with a really sharp cutoff, it's still possible that you will get time-domain artifacts that you don't want. This intuitively makes sense, because you can view the frequency response of the "not-quite-ideal" filter as the response of the ideal filter convolved with a narrow main lobe that smears the response out a bit. Frequency-domain convolution is equivalent to windowing in the time domain; if you look at the impulse response of a very sharp filter that you've designed, it's likely to look much like a $sinc$ function with a window applied to taper the impulse response off at a faster rate than the $sinc$ does on its own.

I give the same advice as with most filter-design problems: you should really try to tackle the problem as quantitatively as possible. Using intuition on what might seem like a good approach can often steer you down the wrong path. Instead, think about the following:

Where is my signal of interest in the spectrum?
What unwanted signals are in the spectrum? Where are they?
How large are the unwanted signals relative to the signal of interest? In order to accomplish my ultimate goal, how much must it be suppressed?
How much distortion can I tolerate on my signal of interest (both in amplitude and phase)?
What computational limitations do I have?

The first four questions will give you filter performance specifications that you can use with any number of design methods to arrive at filters that will achieve your goals. The last question is also important, and can be used to choose between different filter topologies (i.e. FIR versus IIR) to find one that is implementable for your application.


Zero State and Zero Input Responses from Steady State Response

Some books prefer different notations, but I will assume that zero-state response is the one due to pure input $x[n] = \cos(0.2 \pi n) u[n]$ under system zero-initial conditions, while the zero-input response is the one with $x[n] = 0$ and initial state being $y[-1] = 2$.

Now, first the zero-input ($x[n] = 0$ ) response, purely due to the nonzero initial condition ($y[-1] = 2$), of the system can be found by a simple recursion of the causal system as:

$$
\begin{align}
y[0] &= 0.9 \cdot 2 \\
y[1] &= 0.9^2 \cdot 2 \\
y[2] &= 0.9^3 \cdot 2 \\
... &= ... \\
y[n] &= 2 \cdot 0.9^{n+1} u[n] \\
\end{align}
$$

which yields for the zero input response (for $n \geq 0$) as

$$\boxed{ y_{zi}[n] = 2 \cdot (0.9)^{n+1} u[n] }$$

Now we shall find the zero-state response, assuming zero initial conditions. This will be lengthy because of the input form $x[n] = 20 \cos(0.2 \pi n) u[n]$, we can follow the following approach.

Let $X(w)$ be the Fourier transform of the input, and $H(w)$ be the Freuqency response of the system described by the above differential equation, then the output $y[n]$ will have the Fourier transform given by $$Y(\omega) = H(\omega) X(\omega).$$

Now $$H(\omega) = \frac{0.1}{1- 0.9 e^{-j \omega} } $$ and you have to find out what $X(w)$ is, then multiply them to find out $Y(\omega)$ and perform an inverse Fourier transform, after a simple partial fraction expansion. Denoting this solution as $y_{zs}[n]$, then your total solution is

$$y[n] = y_{zi}[n] + y_{zs}[n]$$

Now, try this on your own and see if you can complete the solution... Then write a matlab program to check your result.


As long as the graph consists only of delays, add/subtract and multiply-with-a-constant, then the system will be LTI.

In this case, if the coefficients $k_1 ... k_m$ don't change with time, the system is LTI.


Need a 2 parameter ranking system for image compression

I am not sure if this is going to end up having 2 parameters specifically, but there are a number of things to note:

score depends on two parameters right now, but it is clearly unable to resolve the differences between different schemes. You can add a third parameter which captures the bit planes you are using. This could be either via:

Normalising the number of bit planes used (e.g. $\frac{3}{8} = 0.375$)
but this might still be not enough as a simple reshuffling of the bitplanes can show;
Indexing your score by the number that results from the bitplanes used. So, if you used bitplanes 0,3,5, that would be 2^0 + 2^3 + 2^5, or simply a tuple of $(0,3,5)$ so that implicitly you have the number of bitplanes and which ones were used. This would give you a lookup of (bitplanes_used, score).

Now, the danger you are running here is formulating some $score = f(\Theta)$, where $f$ is a function whose parameter set $\Theta$ keeps growing, more than the 2 you specify. What if we include the ratio of high frequencies to low frequencies in there too? What if we add an entropy related term?...and so on. Ultimately, you might end up with a very large model that is still not informative.

A common thing to do in this case is to penalise large models. This is an extra factor that increases with the number of parameters used. For example, in your case, you could say something like "If two bitplane combinations give me the same score, then I will select the combination that uses fewer bitplanes than the other"

But, it is possible that after all this, your lookups still contain equal scores. In that case, you can collect a subset of the values (that satisfy some criterion in your application) and then randomly select one of them, since all of these values represent a choice of equal score.

Hope this helps


Find Short Clip of Audio Within Longer Clip Of Audio

As good as AUDFPRINT is, I think it is overkill for your problem, as is solving a different problem (where the audio being compared can be arbitrary). You already know what that your shorter file and your longer file capture the same audio scene and you only need to calculate the "offset" between them.

This thread on cross-correlation has a lot of good info:
Cross correlation vs FFT for finding phase between 2 signals

How do I implement cross-correlation to prove two audio files are similar?

If you go the cross-correlation route, I think that phase-correlation (variant) will give you much better results, than simply the basic cross-correlation.

The primary difference is that you normalize your fft coefficients to unit magnitude prior to your inverse-fft, so the phase-correlation is based only on phase information and is insensitive to changes in magnitude/intensity (which makes it more robust, able to synchronize quiet and loud recordings of the same scene)


The covariance matrix of a Kalman filter is a function of the $ Q $ and $ R $ matrices of the model.
If you use a model where $ R $ and $ Q $ are time invariant or known in prior then the calculation of the covariance matrix $ P $ can be done offline and isn't a function of the measurements.

In some cases, advanced implementations of Kalman Filter estimate the covariances $ R $ and $ Q $ on line according to some data gathered in the process of calculating the result.


Phase shift keying is linear modulation.

Digital phase modulation is not necessarily linear.


Can I model process noise as a known "error" in my dynamics while designing a Kalman Filter?

The Extended Kalman filter (EKF) effectively does that: Unmodelled nonlinearities are taken account of by assuming higher process or measurement noise, depending on where the nonlinearity appears.


Delayed impulse does not delay spectral response

The convolution theorem ("Multiplication in time is Convolution in frequency") states:

$$\mathcal{F}\{x\cdot y\}=\mathcal{F}\{x\}*\mathcal{F}\{y\}$$

where $x, y$ are two time-domain signals, $\cdot$ denotes element-wise multiplication and $*$ is the convolution operator. So, the convolution in the frequency domain, is carried out between the Fourier Transforms of the original signals.

As an example, let us use your impulses with frequency response:

$$\mathcal{F}\{\delta(t)\}=1$$

Further, note the time-shifting property of the FT:

$$\mathcal{F}\{x(t-t_0)\}=\mathcal{F}\{x(t)\}\exp(-j2\pi ft_0)$$

where $x(t)$ is any signal you like.
Now, you can readily explain the results you obtain for your shifted impulses:

$$\mathcal{F}\{\delta(t-1)\}=\mathcal{F}\{\delta(t))\}\exp(j2\pi f)$$

where $f=[0/N, 1/N, 2/N, \ldots, (N-1)/N]$ (due to the discrete FT).  Similar for the two-shifted case:

$$\mathcal{F}\{\delta(t-2)\}=\mathcal{F}\{\delta(t)\}\exp(-j2\pi f2)$$

Note that $\mathcal{F}\{\delta(t)\}=1$. You display the absolute value of your spectra: $\left|X(f)\exp(j2\pi ft_0)\right|=|X(f)|$ for any $X(f)$ you like. That's why you see three times the same output. Try plotting real and imaginary part separately, then you will see the effect of time-shifting. You might also have a look here which is an article I wrote about the effect of the time-shift and frequency-shift.


Where can I get an open source implementation of the Parks-McClellan optimal FIR filter design algorithm?

There is an open-source implementation of Parks-McClellan (also known as the Remez exchange algorithm) in GNU Octave, a free-software implementation of a MATLAB-like environment. The function, called "remez", is contained in the "signal" package, which is hosted at Octave-Forge. If you download the package, you'll find "remez.cc", a C++ implementation of the algorithm.

One nice thing about Octave is that it is almost code-compatible with MATLAB, so you can easily port code over to use there if you like. It's a good way to peek under the hood at implementations of algorithms that are provided in MEX form in MATLAB.


The short answer is: Complex exponentials are eigenfunctions of an LTI system and the IFFT modulates different complex exponentials with the QAM data.

More elaborate: The wireless channel can (if it is time-invariant for the duration of one OFDM symbol) be modeled as an LTI system with impulse response $h[n]$ (in discrete-time). So, the received signal is the convolution of the input with $h[n]$. At the receiver, you want to know the data, so if you did directly transmit the data, you would need to perform a deconvolution at the receiver.

One can understand the FFT/IFFT to exploit the convolution theorem (convolution in time is multiplication in frequency). We transmit a signal in time domain at the transmitter, and transform it into the frequency domain at the receiver. Then we know, that the frequency response of the channel has multiplied the frequecy domain of the transmit signal. So, it becomes obvious that we directly define the frequency domain of the transmit signal with our QAM symbols.

This is not the whole story (circular convolution with the impulse response must be achieved such that the convolution theorem holds in discrete, and we use a Cyclic prefix to obtain it).

I've written some articles about OFDM, namely an OFDM walkthrough example and about the CP . Maybe these articles can help you understand better. I've always wanted to write something about the math behind OFDM, but didn't find the time yet.


What is the technical term for the product of "number of pixels" and "number of channel per pixel"?

number of channel per pixel

This is typically referred to as colour depth or bit depth and expressed as bits or 'bits per pixel' (bpp)

  number of pixels

There is no special terminology for this, it's just "Pixel count" or similar. It's unusual to need to refer to the total pixel count of an image without any other context, the dimensions are usually used instead.

  If I calculate 3mn, what is the name of quantity for this value?

There is no special word for this quantity, but the units would typically be bits or bytes, for example

1024 x 768 x 32bpp

= 3,145,728 bytes (3Mb)


Determine the Unit Impulse Response of ARMA Filter for 8 First Indices

As you wrote correctly the function is given by:

$$ y \left[ n \right] = \frac{1}{2} y \left[ n - 1 \right] + \frac{1}{8} y \left[ n - 2 \right] + \frac{1}{2} x \left[ n - 2 \right] $$

Assuming initial conditions of 0 we can easily see that for unit impulse (Namely all zeros but $ x \left[ 0 \right] = 1 $:

$$ \begin{align*}
y \left[ 0 \right] & = 0 \\
y \left[ 1 \right] & = 0 \\
y \left[ 2 \right] & = \frac{1}{2} \\
y \left[ 3 \right] & = 0.5 * y \left[ 2 \right] = \frac{1}{4} \\
y \left[ 4 \right] & = 0.5 * y \left[ 3 \right] + 0.125 * y \left[ 2 \right] = \frac{3}{16} \\
y \left[ 5 \right] & = 0.5 * y \left[ 4 \right] + 0.125 * y \left[ 3 \right] =
\frac{1}{8} \\
y \left[ 6 \right] & = 0.5 * y \left[ 5 \right] + 0.125 * y \left[ 4 \right] = \frac{11}{128} \\
y \left[ 7 \right] & = 0.5 * y \left[ 6 \right] + 0.125 * y \left[ 5 \right] = \frac{15}{256} \\
\end{align*}
$$

Since it has IIR Component it will go on infinitely.

MATLAB Simulation:

The full code can be found in my StackExchange Signal Processing Q43918 GitHub Repository.


Extending 1D window functions to 3D (or higher)

this question is old, but  I needed same function.

I wrote a function that receive one-dimensional window as parameter and transformed into a cube window

function [win] = window3D(w)

    L = size(w,1);
    m1 = w(:)*w(:).';
    win1 = repmat(m1, [1 1 L]);
    m2 = w(:)*ones(1,L);
    win2 = repmat(m2, [1 1 L]);
    win2 = permute(win2, [3, 2, 1]);
    win = win1.*win2;
end


Adjusting corner frequency to constrain maximum disturbance in a high-pass filter

You set $|G/G_0| = 0.001$ and $f = 1kHz$ and then solve your equation (in magnitude form) for $f_c$


So, I myself have taught based on material that uses that scheme, "P/S" and "S/P" after and before the transforms.

Personally, it's nonsense.

What the author tries to say is:

  The IFFT is a mapping of sample vectors to vectors. So, you need a vector as input, not a stream of samples.

What they instead say is:

  We use the terminology from very basic digital logic to deal with complex values.

That's not only a bit awkward, if you asked me, but also inaccurate:

modern FFT implementations in hardware actually do take in stream data, so that operation isn't there in practical hardware implementations
in software implementations, you basically never even deal with samples coming in serially – they always appear en block, in some memory location, so nothing's "serial" to even begin with
S/P nor P/S even say that you should be fully accumulating one vector of $N$ samples, then move on to the next vector, and have zero overlap between these – when applied to FIFOs (which is where the term comes from), that's usually not how they work.

So, simply think of S/P as "get $N$ samples, present these $N$ to the next block as vector".

And P/S is "take this $N$ long vector, and give one sample after the other".

So, these blocks do exactly nothing to your signal - it's just a reinterpretation, if you will, between things that are logically "one after the other", and things that are logically "a vector of $N$ values".


What are the units of my data after an FFT?

In short: if you directly use the voltage output, your FFT amplitude should be in mV as well. If you have the sensor calibration curve, the FFT amplitude should be in Teslas per second (T/s), as you look at a derivative.
If you look at a power density spectrum (squared), then the units above should be squared as well (mV$^2$ or T$^2$/s$^2$).
In general, for a given  unit U, albeit in volt (V), tesla (T), whatever. The FFT sums samples in the original units (U) multiplied by unitless complex values (due to discretization). The units after FFT remain the same as for the original signal, i.e. U.
If you take the absolute value, the same again. For instance, the $0$-frequency index gives you the DC, or average value of your signal  (or at least it is proportional, with a factor of the number of samples, or its square root depending on the normalization).
The continuous Fourier transform, in contrast, "sums" the samples times complex cisoids times a time differential (the $\mathrm{d}t$). Its units would be U.s (s for seconds). Unless you normalize it by something with a unit.


Proper method for generating Channel coefficients in MATLAB

Complex channel coefficient is just a way to represent the independent real coefficients. You just need to generate h = [h_0, h_1, h_2] = [hR_0, hR_1, hR_2] + 1i * [hI_0, hI_1, hI_2]. The independence/correlation between coefficients depend on your model. And if I were not wrong, the number of element of h is the order of your MA model.
The idea behind complex representation is that if the orthogonal set of vector space of baseband signal of bandwidth $[0,W]$ is $\left\lbrace \phi_k(t) \right\rbrace$, the orthogonal set of vector space of passband signal $[f_c-W,f_c+W]$ is $\left\lbrace \phi_k(t) \cos(2\pi f_ct), \phi_k(t) \sin(2\pi f_ct) \right\rbrace$. More details can be found at Read chapter 2.
Furthermore, as real signals have symmetric spectrum, it is convenient to represent these signals by its Baseband Equivalent versions. Thus, a complex representation is used.

So, if you work in the baseband version of signal, channel coefficients are generally complex.

The continuos time white noise is a model of infinite-dimensional signal. When you project this signal to the orthogonal set to obtain the discrete time noise samples, you follow the dimension of the orthogonal set, thus the baseband noise samples are generally complex.

The input data can be complex or real depending on how you want to use available dimensions.

After all, yes, you can process your signal as you have written, regardless real or complex assumptions.

for k = 2:N y(k) = x(k) + 0.6*x(k-1) + 0.3*x(k-2);


I read your question a little differently than the other people who have answered.  The role of the Fourier transform in physics is as a bi-directional one-to-one mapping between conjugate variables.  The most famous example of a conjugate variable pair is position/momentum, but there are many others. Some examples of conjugate variable pairs:

position and momentum over Planck's constant $[x,p/\hbar]$
time and energy $[t,E]$
aperture and the sine of the diffraction angle divided by the wavelength
$\left[x,\frac{\sin\theta}{\lambda}\right]$

One of the consequences in quantum mechanics is that the operators (quantum mechanics uses differential operators which are related to the classical variables) do not commute.  So, using the position/momentum example; $xp\neq px$.  This leads to an uncertainty relation between conjugate variables.

The usual statement of the Heisenberg uncertianty principle is that the more well one knows the position of a quantum particle, the less one can know about the momentum.  Mathematically;
$$
\sigma_x\sigma_ p\geq\frac{\hbar}{2}.
$$
However, the more general statement of the uncertainty principle is that the uncertainty between any pair of conjugate variables is given by
$$
\sigma_1\sigma_2\geq\frac12.
$$
This uncertainty relation leads to some rather strange effects in quantum physics.  For example, the uncertainty relation between energy and time allow for the violation of energy conservation on very brief timescales, a law which is sacred in classical physics.

If you are looking for more information on the uses of the Fourier transform in physics and engineering with many worked examples, have a look at this free textbook by J. F. James.


$2\pi$ periodicity of discrete-time Fourier transform

The argument does not work in continuous time. In discrete time the argument is that
$$e^{j\omega n}=e^{j(\omega+2\pi)n},\qquad n\in\mathbb{Z}\tag{1}$$
This is true because by definition $n$ is an integer. In continuous time we generally have
$$e^{j\omega t}\neq e^{j(\omega+2\pi)t},\qquad t\in\mathbb{R}\tag{2}$$
because $t$ is a real variable. Eq. $(2)$ only holds with equality for integer $t$, but not for all other values of $t$ (which would be necessary for $e^{j\omega t}$ to be $2\pi$-periodic).
But the fundamental reason why the DTFT is periodic is the fact that we deal with sequences, and if we assume that the sequence is the result of sampling a continuous-time signal, then it becomes obvious that different continuous-time signals can result in the same sampled sequence. E.g., take a continuous-time cosine
$$x(t)=\cos(2\pi f_0 t)\tag{3}$$
Now we sample $x(t)$ at instances $nT=n/f_s$, where $f_s$ is the sampling frequency:
$$x[n]=\cos(2\pi f_0 nT)=\cos(2\pi n f_0/f_s)\tag{4}$$
I leave it as an exercise to you to show that the frequencies $f_0=f_s-\Delta f$ and $f_0=f_s+\Delta f$ (with some $\Delta f>0$) result in the same sequence $x[n]$. This phenomenon is called aliasing. This means that frequencies above $f_s/2$ (also called the Nyquist frequency) cannot be distinguished from the corresponding frequencies below the Nyquist frequency. That's why the DTFT is periodic. There cannot be any new spectral information outside the frequency interval $[-f_s/2,f_s/2]$. If you haven't learned it yet, check out the sampling theorem, which will also help you understand the nature of spectra of discrete-time signals.
Also have a look at this answer to a related question.


What preprocessing technique can help me to filter the background text noises using matlab...the image is shown below?

This problem is called (ink) bleeding. A popular approach is to separate or segment the document image using graphical models. Relevant papers include:

Document Ink Bleed-Through Removal with Two Hidden Markov Random Fields and a Single Observation Field
User-Assisted Ink-Bleed Reduction
Bleed-through removal in degraded documents
Directed Assistance for Ink-Bleed Reduction in Old Documents
Ink-Bleed Reduction Using Functional Minimization


Interpolating irregularly missing data points of regularly spaced data

I found a method similar to the zero-padding FFT interpolation but in which the missing samples can be nonuniformly spaced. I think it exactly solves the problem you have. The method is described here:

J. Selva, "FFT interpolation from nonuniform samples lying in a regular grid", IEEE Trans. on Signal Processing, vol. 63, n. 11, June 2015.

The code is available here.


What is the best way to show beamformer enhances signal of interest over interference

There's an article that defines some  measures for blind source separation that I used to measure the output of a beamformer Performance Measurement in Blind Audio Source Separation. Vincent E., Gribonval R. and Févotte C. Transactions on Audio, Speech, and Language Processing, Vol. 14, No. 4, July 2006. They basically use orthogonal projections to decompose the output signal and use these projections to calculate Signal to Noise, Interference, Distorsion and Artifacts ratio, I know the paper is based on BSS, but it could be an starting point.

http://www.irisa.fr/metiss/demonstrations-fr/bssperf.html


By estimating F1 as closely as you can and removing it from your signal.  Then analyze the residue for F2.

This is how I would do it:

1) Partition your signal into small windows that are about 2 1/2 or 3 1/2 cycles of F1 in length.  Since you frequency is varying, the window sizes will too.

2) Apply this technique in each window to get the paramters for that window.

A Two Bin Solution

3) Construct your best fit of F1 by constructing a signal using the parameters for each window.  Extend the construction beyond your windows and do a sliding linear scale on the overlaps to merge them.  This is the more efficient approach.

3) (Alternative)  Derive an unevenly samples sequence from the parameters, curve fit them to make signals with varying frequency.  Several ways to do this, all require more calculations, and probably aren't worth it in a noisy signal.  If you have high SNR and looking for high precision, this is better way.

4) Subtract the constructed signal from your original signal.

5) Do a FFT on the leftover on convenient window sizes.

5) (Alternative) If you can now recognize (or have an estimate of) F2's frequency), apply the same technique as in step 2.  Calculating two DFT bins is a lot fewer calculations than doing a full FFT.

If F1 is varying slowly, you may want to use a larger number of cycles + 1/2.  The longer the window the better the noise mitigation, but also the worse fit if it is varying too much.

On second thought.  It would actually be easier to pick windows in step 1 to be whole cycles of F1 in length, then take the FFT.  This will minimize the leakage from F1 and you can analyze F2 independently since it is many bins away.  The same formulas for "The Two Bin Solution" will work for that.  If F2 happens to be really close to an integer number of cycles per frame  (within 0.0000001% or so), use the single bin value directly instead or adjust your window length by including a few more F1 cycles to make F2 near a whole plus a half number of cycles.  The two bin formulas are the most robust (noise resistant)  when the frequency is halfway between the bins.  This will give you a better read than a single bin on whole cycles.


expected value of two LTI output signals multiplied (cross correlation)

Assumimg that the two linear systems are BIBO-stable, the random processes $\{Y_1(t)\}$ and $\{Y_2(t)\}$ are zero-mean WSS Gaussian processes with autocorrelation functions and power spectral densities given by
\begin{align}
R_{Y_1} &= \sigma^2 (h_1 \star \tilde{h}_1)\\
R_{Y_2} &= \sigma^2 (h_2 \star \tilde{h}_2)\\
S_{Y_1} &= \sigma^2 |H_1|^2\\
S_{Y_2} &= \sigma^2 |H_2|^2
\end{align}
In fact, the processes are also jointly Gaussian and jointly WSS processes with cross-correlation function $$R_{Y_1, Y_2}(\tau) = E[Y_1(t), Y_2(t+\tau)] = \sigma^2 (h_1 \star \tilde{h}_2)$$ and cross-power spectral density
$$S_{Y_1,Y_2}(f) = \sigma^2 H_1(f)H_2^*(f).$$
The OP wants to find $E[Y_1(t)Y_2(t)]$ which is given by
\begin{align}E[Y_1(t)Y_2(t)] &= R_{Y_1,Y_2}(0)\\
&= \sigma^2 h_1\star \tilde{h}_2\big|_{\tau=0}\tag{1}\\
&= \sigma^2\int_{-\infty}^{\infty} H_1(f)H_2^*(f) \,\mathrm df\tag{2}
\end{align}
since the OP prefers the frequency-domain calculation. Personally, given only $h_1$ and $h_2$ (and not $H_1$ and $H_2$), I would say that it is easier/cheaper to just grind out $(1)$ rather than use the frequency-domain calculation $(2)$ but the OP might have specific reasons for opting for the frequency-domain calculation. In particular, $\operatorname{var}(Y_1+Y_2)$ which is what the OP seems to really want to find is just
$$\operatorname{var}(Y_1+Y_2) = \sigma^2\begin{cases}\displaystyle\int_{-\infty}^\infty |h_1(t)+h_2(t)|^2 \,\mathrm dt,\\
\displaystyle\sum_{n=-\infty}^\infty |h_1[n]+h_2[n]|^2\end{cases}$$
which seems easier than the frequency-domain version of the same calculation, but ymmv.


If we define samples manually how could there be a Sampling Frequency?

Units can be arbitrary, such as ISO Meters, or the length of the current king's foot (maybe longer than last year after the assassination of some shorter king).

A good default unit if there isn't some ISO standard, useful natural unit, or you can't think of one to make up (how long is your monitor's display, or your left foot, etc.), is the length of your window of data.  e.g. 1 cycle per your vector of 256, 2 cycles, 3 cycles, etc. per your vector length.  If your vector displays as 3 inches long on your display, then you can also easily convert cycles/vec into cycles/inch, or per/mm, per/fortnight, etc.


Detecting and removing interferences from a signal

I'd recommend using a median filter to smooth your signal. A median filter will get rid of the outliers, or the spikes in your signal. The length of the median filter will have to be determined by how frequently the spikes appear in your signal. The median filter is available in Python as numpy.median.
If you want to retrieve the spikes (the interference signal), you can simply subtract the filtered signal from the noisy signal to get the interference.


i am a little confused by the first sentence. i dunno how a wavetable oscillator scans a sound file at all. some algorithm can scan a file and extract either grains (or little "wavelets", which is a term i have misused in the past) or single-cycle wavetable "snapshots". but the wavetable oscillator starts with wavetable(s) and outputs a periodic or quasiperiodic tone.

about aliasing, if you extract wavetables from a real tone (and let's assume you have a decent pitch detector, so you know how long each period is to a precision of fraction of sample, and decent interpolation, so you you can resample the period to a, say, 256 or 1024 point wavetable) and you up-pitch the tone by more than a couple of semitones, there is a danger of aliasing unless you remove some of the higher harmonics. you can remove these higher harmonics by performing a DFT (a.k.a. FFT) on the wavetable and zero some of the bins centered around bin N/2 . then inverse DFT back to a wavetable.

actually, the number of bins centered around bin N/2 that are already zero is dependent on the degree of upsampling in the interpolation when the period in the soundfile is converted into a wavetable.  consider this example:  suppose the sample rate of the sound file is 44100 Hz and you have a note that is at middle C, with a fundamental at 261.6255653 Hz.  this means that your period length is 168.5615 samples and a really good pitch detector can return a period value with that precision (even though i doubt it would be that accurate, precision and accuracy are different things).  now a 128-sample wavetable would be insufficient to represent that waveform, or at least possibly insufficient to represent all of the harmonics of that waveform.  there are conceptually 84 harmonics.  a 256-samplewavetable would be sufficient to represent all of the harmonics, but let's say, for playback reasons, you resample that 168.5615 sample period to a 1024-sample wavetable, $x[n]$.  that means you have a step size (the reciprocal of the sample-rate-conversion ratio) of 168.5615/1024 = 0.16461084 .

now, if you were to FFT that 1024-sample wavetable, you would possibly have non-zero bins for $X[0]$ up to $X[84]$ (and also the mirror image of $X[1024-84]$ up to $X[1023]$ .  all of the bins from $X[85]$ to $X[1024-85]$ should be zero if the interpolation kernel was any good.

now, when you play that wavetable back at the same note, middle C, then the sample stride would be the reciprocal of the step size above, or 6.074934 wavetable samples per output sample (assuming the same 44100 Hz sample rate).  and you would have no aliasing because all of those harmonics above the 84th harmonic are zero.  but if you start pitching that note higher, your stride will be longer.  if you pitched it 7 semitones higher (to G), then the stride would be 9.1021167 wavetable samples per output sample.  but then your 84th harmonic gets bumped up to around 33 kHz which is above Nyquist and you will have aliasing.

so what you need to do, if you plan to pitch that note up 7 semitones (i would not pitch a sampled piano up that far, because it will sound "plinky"), you must lose all of those harmonics above the 56th harmonic.  ($56 \approx 84 \times 2^{-7/12} $).  you do that with the FFT, zeroing bins, and then inverse FFT


There is a difference between a quantized signal, i.e. taking only integer values, and a sampled signal, i.e. the value of which is only looked at on spaced instants.

A physical signal is usually continuous in values and in time.

An A/D converter turns a physical signal in a quantized one; if the converter is asynchronous, it generates a time-continuous signal.

Discrete-time/continuous-value signals are created by sample & hold circuits.

And digital signals processed by computers are usually both quantized and time-discrete.


identify equilibrium region after large transient

It is difficult to think about robustness without knowing anything about the noise characteristics.

For example, it might be enough to obtain a number of these waveforms (20-50), obtain a box plot of the amplitude at each time instance and use that information to derive a very simple threshold below which the signal enters a region that can be safely characterised as equilibrium.

If there is enough noise to cause the signal to "jump" above the threshold, then you might want to try introducing some form of hysteresis in the Thresholding process. In other words, you might consider "equilibrium" the point when the waveform has dropped below a certain threshold and has stayed there for at least a $u$ amount of time.

If the noise characteristics are such that it is impossible to "tune" such a threshold then a simple 3-5 tap one-dimensional median filter or a very simple low pass filter might be enough, to smooth out the waveform sufficiently and enable you to apply a threshold more reliably.

If the process is such that the waveform is "floating" and a fixed threshold is impossible, then you would be looking at the point when the rate of decay is sufficiently flat. Since you already know that you have an exponentially decaying waveform, you should use this to your advantage. In this case, you can still apply some sort of smoothing pre-processing (e.g. filtering, if required) and then fit the data to an exponentially decaying model. You can then obtain the first order finite difference and set a threshold on that variable for when it drops sufficiently close to zero denoting that the curve is flattening out.

But in your case, you seem to have some sort of undershoot as well (some sort of inertia effect (?)). In extreme cases, the exponential decay model might not be a good assumption in terms of fitting the exact waveform. In that case, it might be better to do a DFT on the settling segment and try to quantify the settling time by known results from automatic control system theory.

Hope this helps.


How to compute beam pattern given a delay-and-sum filter in spectral domain?

One of the method for calculating beam pattern could be, you simulate/capture audio from different directions on array and filter these captures with filters($H$'s) generated. the audio coming from direction of beam pattern will be passed as it is, but the audio coming from other directions will be attenuated. the amount of attenuation with respect to spatial location is nothing  but the beam pattern .
for implementation,
Take audio in time domain, window it, take n-point DFT of this windowed time domain audio after this do element wise multiplication with filter coefficients to get beamformed output.


You are mostly correct.  The equation should be $w(t) = sin(2\pi\frac{1}{20}t)$.  The variable in the sin should be $t$, not $f$.  The frequency is determined by the $\frac{1}{20}$, not by the variable.

Other than that you are correct.  The period is 20s, which makes it a $\frac{1}{20}$ Hz sin wave.


DFT have independent variable in between 0 to $\frac{1}{2}$ no of samples

This statement is only true if your time-domain signal $x[n]$ is real. A real signal has a Fourier transform that is conjugate-symmetric. Stated differently the "negative frequencies" are just the conjugate of the corresponding "positive frequencies." Therefore, the negative frequencies add no new information, so they can be discarded from the DFT of a real signal.


Extracting audio track to determine the song beat

Sonic Annotator. Good: their beat tracking algorithm is the most sophisticated. Bad: this is not easy to integrate and interface into your own application (you need to load their own flavor of plug-ins called "VAMP"); and the beat tracker is closed source - available as a pre-compiled binary lib.

Aubio. Good: very easy to integrate, plain C with few dependencies. Bad: the quality of the beat-tracking is so-so.

Marsyas. Good: beat-tracking accuracy is OK (but not as good as Sonic Annotator's). Bad: many dependencies, very complex framework.

Depending on what is your target application (is the audio to analyze a recording of a well-known/distributed song, and does your system have internet access), you could also use the Echo Nest APIs.

If you want to roll your own, I suggest you this as a starting point - start porting this rather simple matlab implementation to your language of choice, and you'll be at 80% or 90% of the accuracy of the most sophisticated solution.


IIR coefficients that simply boost the signal

I will assume that your IIR filter is an order-2 biquad. But the reasoning remains the same if the implementation is different.

Original filter : $H(z) = \frac{b_0+b_1\cdot z^{-1}+b\cdot z^{-2}}{1 + a_1\cdot z^{-1} + a_2\cdot z^{-2}}$

If you want to boost your signal by G

$H_G(z) = G\frac{b_0+b_1\cdot z^{-1}+b_2\cdot z^{-2}}{1 + a_1\cdot z^{-1} + a_2\cdot z^{-2}}$

simply distribute G to the numerator

$H_G(z) = \frac{G\cdot b_0+G\cdot b_1\cdot z^{-1}+G\cdot b_2\cdot z^{-2}}{1 + a_1\cdot z^{-1} + a_2\cdot z^{-2}}$

Simply multiply your b0-b2 coefficients by G.

In the case where you simply want to boost the signal without filering it at all. Set $b_0$ to $G$, $b_1$, $b_2$ to $0$, as well as  $a_1$ and $a_2$ to 0.


How to delay an audio signal by 1 microsecond?

It's called a "fractional delay".

A really comprehensive analysis on how to do this can be found here: https://ieeexplore.ieee.org/document/482137/ and the design tools are still available here http://legacy.spa.aalto.fi/software/fdtools/


What is impulse response in simple AWGN channel?

A linear system cannot add anything new to its input signal, so additive noise is modeled separately from any linear distortion of the signal. In your case with just scaling and delay, the channel impulse response is $h(t)=k\delta(t-\Delta t)$, and the complete channel is modeled by first filtering the input signal with an LTI system with impulse response $h(t)$ and then adding the noise.


How determine the delay in my signal practically

The Least Mean Square solution to find the "channel" or response of the filter is provided by the following MATLAB/Octave Code using  the input to the filter as tx and the output of the filter as rx. For more details on how this works, see this post: Compensating Loudspeaker frequency response in an audio signal:

function coeff = channel(tx,rx,ntaps)
    % Determines channel coefficients using the Wiener-Hopf equations (LMS Solution)
    % TX = Transmitted (channel input) waveform, row vector, length must be >> ntaps
    % RX = Received (ch output) waveform, row vector, length must be >> ntaps
    % NTAPS = Number of taps for channel coefficients
    % Dan Boschen 1/13/2020

    tx= tx(:)';   % force row vector
    rx= rx(:)';   % force row vector
    depth = min(length(rx),length(tx));
    A=convmtx(rx(1:depth).',ntaps);
    R=A'*A;       % autocorrelation matrix
    X=[tx(1:depth) zeros(1,ntaps-1)].';
    ro=A'*X;      % cross correlation vector
    coeff=(inv(R)*ro);   %solution
end

The case the OP uses of a comb filter is one of the most challenging as it depends on signal energy at each frequency for a solution (this is specifically why linear equalizers, which this function is doing if you swap rx and tx, do not perform well in frequency selective channels and result in noise enhancement at the null locations).  Below the frequency response of the test filter as determined with MATLAB or Octave showing the multiple frequency nulls associated with such a comb filter:

b=[1,zeros(1,39),-1];
freqz(b,1,2^14)     % 2^14 samples to show nulls

MATLAB or Octave script to demonstrate use of above function and determine the delay between output and input:

%% Filter with OP's example
b=[1,zeros(1,39),-1];     % numerator coefficients
a = 1;                    % denominator coefficients

%% Generate signal using OP's code
Fs=1000;
t=0:1/Fs:(5-1/Fs);
wi=blackman(length(t))';
rn=+rand(1,length(t))*.2;
x1=sin(2*pi*13*t).*wi +rn;
x2=sin(2*pi*25*t).*wi +rn;
x3=sin(2*pi*2*t).*wi +rn;
x=[x1,x2,x3];

% Filter
y=filter(b,a,x);

%% Test filter estimation
cf=channel(x,y,61);

%compare original and estimated channel
subplot(2,1,1)
stem(b)
title("Actual Channel Response")
xlabel("Sample Number")
subplot(2,1,2)

stem(cf)
title("Estimated Channel Response")
xlabel("Sample Number")

We could have used exactly 41 taps when calling the function channel and it would have resolved the solution properly; however it is best practice to start with a much longer filter length, evaluate the result and then reduce the taps accordingly. In actual practice under noise conditions using more taps than necessary will lead to noise enhancement, so it is good to do the final solution with the minimum taps necessary to capture the dominant tap weights.

Observe with the groupdelay plot using MATLAB and Octaves grpdelay command the issue with not being able to resolve the delay where no signal gets through the filter (it would be difficult for anything to determine the delay of a single tone at one of these frequencies that is nulled by the filter!) but it is able to accurately determined the delay where signal energy exists. Similarly the waveform itself must have energy at all frequencies where we are seeking a solution. The spectral density of the OP's test waveform was sufficiently spread out over all frequencies to be suitable for this purpose (this is the reason why pseudo-random waveform make for good "channel sounding" patterns).

This plot is for comparison to the OP's plot showing that the group delay for this filter is 20 samples.

Here is another test case showing how this works with a more reasonable channel (no deep nulls) with the following coefficients and frequency response below that:

b = [0.2 .4 -.3 .4 .3 .1];

The solution is indistinguishable between actual and estimate, so I added noise to make it more interesting, using same x and y from code above:

% add noise
noise = 0.351*randn(1,length(y));
yn = y + noise;
snr = 20*log10(std(yn)./std(noise));
%% Test filter estimation
cf=channel(x,yn,10);


Imagine for one second, that you just plotted your daubechies-4 wavelet, as you can
see here in red.

Now imagine that you take this waveform in red, and simply do a cross-correlation with your signal. You plot that result. This will be the first row of your plot. This is scale-1. Next, you dilate your Daubechies-4 wavelet, (that is, you simply make it 'stretch' it in time, by some factor). Then, you again do a cross-correlation of this new waveform with your signal. You then get row two of your plot. This is scale-2. You keep doing this for all scales, which means you keep taking your original 'mother' wavelet, and you keep dilating, then cross-correlating, dilating, and cross-correlating, etc, and you just plot the results one on top of the others.

This is what the CWT plot is showing you. The results of performing a cross-correlation of your signal with a wavelet at different scales, that is, at different dilation (stretch) factors.

So let us interpret your image. In the first row, you can see that you have weak amplitudes in your cross-correlation. That means that it is telling you, almost nothing in your signal correlates, (or 'matches') your wavelet, when it is at scale-1, (when it is at the smallest scale). You keep stretching your wavelet and correlating, and still it is not matching anything in your signal, until you reach say, scale-31. So when you stretch your wavelet 31 times, and performed a cross-correlation, you start to see some bright spots, which means you are getting good cross-correlation scores between your stretched wavelet, and your signal.

If you look however at the very top, we have the brightest spots. So for scale-46, you made that row by stretching your original wavelet 46 times, and then cross-correlated it with your signal, and then that is your row-46. So you see a lot of nice bright spots. You can see that at positions (x-axis) ~25, ~190, and ~610, I have bright spots. So that is telling you, you have some feature in your signal, that very closely matches your wavelet that is stretched 46 times. So you have 'something' at those locations that closely match your wavelet at this scale.

(Of course, in your case, you have used noise, so the positions that I have talked about are random - that is, there is nothing really 'interesting' going on. Do a CWT with a sine pulse and what I am saying can be made clearer to you.)

In summary, the CWT is simply showing you all possible correlation scores between your template/matched filter (in this case daub-4 wavelet), at different positions, (x-axis) , also at different stretched factors, (y-axis).

Hope this helped.


I don't fully understand what the OP means by DQPSK but suspect that he believes that the D in DQPSK stands for "differentially coherent demodulation of" QPSK, that is, DQPSK means differentially-coherent demodulation of a plain vanilla QPSK signal, and that in his terminology, differentially encoded DQPSK is differentially-coherent demodulation of differentially-encoded QPSK. Most people would take DQPSK to mean the latter because the former interpretation of DQPSK (differentially-coherent demodulation of a plain vanilla QPSK signal) is a really bad idea.
Looking only at the output of a QPSK transmitter, one sees an RF signal whose RF phase has constant value over successive $T$-second epochs, but can change abruptly at the transition from one epoch to the next from one value in $\left\{\frac{\pi}{4},\frac{3\pi}{4}\frac{5\pi}{4}\frac{7\pi}{4}\right\}$ to another value in $\left\{\frac{\pi}{4},\frac{3\pi}{4}\frac{5\pi}{4}\frac{7\pi}{4}\right\}$ (or stay the same, too). For some details, see this answer of mine. Assuming that the signal is a plain vanilla QPSK signal, the (coherent) receiver needs to work hard (using PLLs or Costas loops) to acquire both the carrier frequency and the carrier phase of the incoming signal so that the receiver's local oscillator can be synchronized with the incoming signal and coherent demodulation is possible. With such synchronization, the QPSK system essentially works as two parallel BPSK schemes operating on phase-orthogonal carriers, and the I and Q branches of the receiver make decisions independently as to whether their separate matched-filter outputs correspond to a transmitted $1$ or a transmitted $0$. Thus, for a given value of $E_b/N_0$ (the "bit energy to noise density" or BEND ratio as I like to call it), the BER for coherent demodulation of QPSK is the same as the BER for BPSK. Note however, that the QPSK signal is transmitting $2$ bits in $T$ seconds and so the signal power is doubled.
Next, let us consider  differentially coherent demodulation of the plain vanilla QPSK signal. Here it is not necessary to determine the RF carrier phase -- synchronization in just the carrier frequency suffices --  but all that can be determined is whether the RF phase of the current epoch is the same as the RF phase during the previous epoch, or that the current RF phase differs from the previous RF phase by $\frac{\pi}{2}$ or by $\pi$ or by $\frac{3\pi}{2}$. Now, if we know what the previous dibit was ($00, 10, 11, 10$ as the case may be), we could deduce what the current dibit is. But, this notion is fraught with the possibility of not just error in deciding what the current dibit is, but error propagation. If we get the current dibit wrong, the next dibit will be wrong too because its value will be determined from the wrong value of the current dibit, and the next dibit decision will be wrong, and the next ad infinitum until another dibit error of just the right kind fortuitously restores everything back to normal. Moral: Don't try to do differentially-coherent demodulation of plain vanilla QPSK; it doesn't work well at all regardless of what the OP believes.
If one is doing differentially-coherent demodulation of differentially-encoded QPSK, then the information is being carried not in the dibit values during the epochs, but in the change in dibit value between epochs. For example, the data bits for the current epoch are determined by the change in phase that occurred during the transition between the previous epoch and the current epoch. Such a change of phase is precisely what the differentially-coherent QPSK receiver is designed to detect (cf. the previous paragraph). It is worth emphasizing that differential encoding in QPSK is not the same as doing differential encoding of the binary data streams of the I and Q branches; it is different. For details of thhe differential encoding, see this answer. For more information ow soft outputs from the matched filters are processed in a DQPSK receiver, see here.


The resulting spectrum appears to be your Doppler shift (possibly as a ratio to half the sampling rate, depending on what the horizontal axis on your plot represents). Without seeing all the processing that was done to generate the "Filtered Chirp Sequence" I cannot be certain but I assume it is a complex result from the subsequent plots, which would rotate in phase at the rate of the Doppler offset. To confirm this, run the test or simulation with known 0 Doppler or Clock offset and confirm that the resulting phase is constant from sample to sample.


Is there a function in python which returns the Amplitude/Sound Pressure Level of a sound file(like .wav) file?

No, you need to know the full Transfer function from the acoustic response at the sensors face to your A/D converter.  Transducer responses vary with frequency. You need a calibrated reference at a minimum.

A reference sensor will have a response curve from the manufacturer.  An average would be useful as well.


What Resources Are Recommended for an Introduction to Signal Processing (DSP)?

My recommendation in terms of text books is Rick Lyons's Understanding DSP.  My review of the latest edition is here.

I, and many others from the ${\tt comp.dsp}$ community and elsewhere, have helped Rick revise parts of the text since the first edition.

For self-study, I know of no better book.

As an on-line, free resource, I recommend Steve Smith's book.  Personally, I prefer Rick's style, but Steve's book as the advantage of online accessibility (and the online version is free!).

Edit:

Rick sent me some feedback that I thought I'd share here:

For your colleagues that have a copy of my DSP book,
I'll be happy to send them the errata for my book.  All they have
to do is send me an E-mail telling me (1) The Edition Number,
and (2) the Printing Number of their copy of the book.  The Printing
Number can be found on the page just before the 'Dedication' page.
My E-mail address is:
R.Lyons [at] ieee.org

I recommend that your colleagues have a look at:
http://www.redcedar.com/learndsp.htm

Rick also gave me a long list of online DSP references.  There are way too many to put here. I will see about setting up a GoogleDocs version and re-post here later.


Calculating amplitude, time-shift and phase of third sinusoid as a sum of two others

Here's how I would do it.

Use the Cosine angle addition formula on $x_1$ and $x_2$:

$$ \cos( \alpha + \beta) = \cos( \alpha )  \cos(\beta)-\sin( \alpha )  \sin(\beta) $$

Like this:

$$ x_1(t) = a_1 \cos( \omega t + \phi_1 ) $$

$$ x_1 = a_1 \cos( \omega t ) \cos(\phi_1)- a_1 \sin( \omega t ) \sin(\phi_1) $$

The $\omega$s will be the same, so:

$$ x_2 = a_2 \cos( \omega t ) \cos(\phi_2)- a_2 \sin( \omega t ) \sin(\phi_2) $$

You can now rearrange and add them together:

$$
\begin{align}
 x_3 &= x_1 + x_2 \\
 &= \left[a_1 \cos(\phi_1) + a_2 \cos(\phi_2) \right] \cos( \omega t )
- \left[a_1 \sin(\phi_1) + a_2 \sin(\phi_2) \right] \sin( \omega t )
\end{align}
$$

Since you know the $\phi$s and $a$s, you can put this in the form:

$$ x_3 =  C \cos( \omega t ) - D \sin( \omega t ) $$

Now, you just have to reverse the process to get it back into time phase form.

$$ \theta = \operatorname{atan2}(D,C) $$

$$ M = \sqrt{ C^2 + D^2 } $$

Therefore:

$$ x_3 =  M \cos( \omega t + \theta ) $$

$$ x_3 =  M \cos\left( \omega \left( t + \frac{\theta}{\omega} \right)  \right) $$

Now, just plug and chug and you should have it.

The latter part might be better understood in a different order:

$$
\begin{align}
x_3(t) &= a_3 \cos( \omega t + \phi_3 )\\
&= a_3 \cos( \omega t ) \cos(\phi_3)- a_3 \sin( \omega t ) \sin(\phi_3)
\end{align}
$$

From the equations above:

$$
\begin{align}
C &= a_3 \cos(\phi_3)\\
D &= a_3 \sin(\phi_3)
\end{align}
$$

From there it follows:

$$ a_3 = \sqrt{C^2+D^2}$$

and

$$ \phi_3 = \tan^{-1}\left(\frac{D}{C}\right) $$

The thing is, both of these last two equations have two possible solutions each and they need to be matched.  Conveniently, the atan2 function is available on most platforms which yields the correct angle from the latter for the positive root of the former.

This solution is actually the proof of a very important principle:  When two pure tones of the same frequency are added together the result is a pure tone of the same frequency although it is possible for it to have zero amplitude.

In order for the zero amplitude case to occur, the two tones have to have the same amplitude and be a half cycle ($\pi$ radians, 180 degrees) out of phase.  This is called complete destructive interference.

If the two tones are in phase, the amplitudes are merely added and the phase remains the same.  This is technically called complete constructive interference.

By repeated application of this principle, the same holds true for any linear combination of tones.

This shows up in a surprising number of places.


which of the analog sinusoidal frequency can not pass through the filter?

The impulse response of the given transfer function $H(z)= 1 + z^{-1} + z^{-2}$ is $h[n] = \delta[n] + \delta[n-1] + \delta[n-2]$.

This is a rectangular window with three taps and its Frequency response magnitude is given by:

$$ H(w) = \frac{\sin( 3 w/2 )}{\sin(w/2)} $$

It'll have zeros (nulls) at the frequencies given by:

$$ \sin(3w/2) = 0 \longrightarrow 3 w_k / 2 = \pi k \longrightarrow w_k = 2\pi k / 3 ~~~~, \text{for} ~~~k = 0,1,2 $$

Note that the zero for $k=0$ is cancelled by the denominator and therefore there exists only two zeros $w_1 = 2\pi/3$ and $w_2 = 4\pi/3$

Now note that those two zeros are at mirror locations about $w=\pi$ which means that they will be the same, indistinguishable, discrete time frequencies.

Also you shall consider only the first zero , $w_1 = 2\pi/3$ to be the frequency of the analog sinusoidal that will be nulled by the filter.

The given sampling frequency of $f_s = 18$ kHz maps $w_1= 2\pi/3$ to the analog frequency of

$$ f_1 = \frac{f_s w_1}{2\pi}  = \frac{18k 2\pi/3}{2\pi} = 6 kHz$$

The answer therefore is $f_1 = 6$ kHz which cannot pass through your filter when sampled properly at the sampling rate of $f_s = 18 $ kHz. This assumes that there is no aliasing during the sampling process which is prevented by an analogue anti-alasing filter before sampling. If this filter is omitted and sampling is performed then a sinusoid at the frequency of $f_2 = 12$ kHz will be aliased into a frequency of $f_1 = 6$ kHz and therefore it won't also pass through the filter. Note that there will be infinetely many new frequencies which will map into $f_1 = 6$ kHz and won't pass through the filter in such a case.


The block diagram you have posed is composed of serial cascade of three IIR sections each having its own input output relationship. The first two are second order stages and the last one is first order.

If you know about Z-transforms, you could easily find the algebraic relation between those signal variables in terms of their Z-transforms.

I would denote the intermediate summation nodes as $w[n]$ and $v[n]$ together with the input $x[n]$ and output $y[n]$ the Z-domain relations become:

$W(z) (1 - a_1 z^{-1} - a_2 z^{-2}) = X(z)$
$V(z) (1 - b_1 z^{-1} - b_2 z^{-2}) = W(z)$
$Y(z) (1 - c_1 z^{-1} ) = V(z)$

where I have denoted the system coefficients as $a,b,c$ respectively. (Do not mix this $a,b$ with the standard LCCDE expression coefficients $a,b$.)

Then simple algebra tells us the relation between $X(z)$ and $Y(z)$ as:

$$Y(z) = 2 \left( \frac{1}{1 - c_1 z^{-1}} ~ \frac{1}{1 - b_1 z^{-1} - b_2 z^{-2}} ~\frac{1}{1 - a_1 z^{-1} - a_2 z^{-2}} \right) X(z)$$

By first expanding the product polynomials and then mapping back from Z-domain to time domain you can express $y[n]$ as a function of $x[n]$ (and feedback)


Spectrogram of a single tone complex signal has two dark lines?

Since this is a constant spectrogram, you could just as well have just averaged the |FFT|² and plotted that! (The most colorful way of visualizing things isn't always the optimal one; your signal doesn't change over time, so you don't need the time axis of the spectrogram at all.)
Quite possibly, in that "easier" representation, you would have spotted this:
Your spectrogram does FFTs of length 512. That's not a multiple of the period of your discrete complex sinusoid.
Therefore you'll see leakage. (You already know that effect!)
However, to reduce that resolution-reducing leakage effect, specgram applies a Hanning window by default to your 512-sample chunks that you FFT.
That hanning window leads to zeros:
plt.plot(20*np.log10(np.abs(np.fft.fft(fc1[:512]*np.hanning(512)))))

And there you go!


Two successive symbols in the demodulator are $Z_1 = (X_1,Y_1)$ and $Z_2 =(X_2,Y_2)$
where $X$ is
the output of the I branch and $Y$ the output of the Q branch of the receiver. The
hard-decision DBPSK decision device considers the question:

  Is the new symbol $Z_2$ closer to the old symbol $Z_1$ or to the
  negative $-Z_1$ of the old symbol?

and thus compares

$$(X_2-X_1)^2 + (Y_2-Y_1)^2 \gtrless (X_2+X_1)^2 + (Y_2+Y_1)^2$$

which can be simplified to a sign comparison on $\langle Z_1,Z_2\rangle = X_1X_2+Y_1Y_2$.
Note that this is essentially asking

  Are the two vectors $Z_1$ and $Z_2$
  are pointing in roughly the same direction (in which case the inner product
  or dot product is positive) or in roughly opposite direction (in which case
  the dot product is negative)?

A third viewpoint thinks of $Z_1$ and $Z_2$ as
complex numbers and asks

  Is $\text{Re}(Z_1Z_2^*) = X_1X_2+Y_1Y_2$ positive or negative?

The soft decision decision device simply
passes on the exact value of the dot product to the soft decision decoder
which may opt to quantize dot products that are very large in magnitude
into hard decisions and continue waffling on the rest.  This is what
the decision rule stated in the OP's question is, where large is taken as
exceeding $1$ in magnitude.

In DQPSK, the encoding uses one of two conventions:

the signal phase is delayed by $0, \pi/2, \pi, 3\pi/2$ according
as the dibit to be transmitted is $00, 01, 11, 10$
the signal phase is advanced by $0, \pi/2, \pi, 3\pi/2$ according
as the dibit to be transmitted is $00, 01, 11, 10$

Note that a DQPSK signal is not the sum of two DBPSK signals
modulated on phase-orthogonal carriers, but the I and Q bits
jointly affect the net carrier phase.

For demodulating a DQPSK signal, the decision device needs to ask

  Which of the four symbols $Z_1,\quad jZ_1 = (-Y_1,X_1),\quad -Z_1,\quad -jZ_1 = (Y_1,-X_1)$
  is $Z_2$ closest to?

Thus, in addition to the comparison

$$(X_2-X_1)^2 + (Y_2-Y_1)^2 \gtrless (X_2+X_1)^2 + (Y_2+Y_1)^2$$

it is necessary to compare

$$(X_2+Y_1)^2 + (Y_2-X_1)^2 \gtrless (X_2-Y_1)^2 + (Y_2+X_1)^2$$

which works out to looking at $\text{Im}(Z_1Z_2^*)$ in
addition to $\text{Re}(Z_1Z_2^*)$ and making the decision
according as which quantity has the largest magnitude
and the sign of the largest magnitude.  The details
of how the soft-decision decoder uses the decision
statistic $Z_1Z_2^* = (\text{Re}(Z_1Z_2^*), \text{Im}(Z_1Z_2^*))$
will determine how these numbers are further
massaged.


Frequency Modulation to arbitrary Audio signals

I have, in the past, used delay modulation on arbitrary audio signals.  but the modulation waveform itself was not arbitrary, but was a sinusoid.  but we could speed this modulation oscillator up to an audio rate.

Now that I think of it, I tried hooking that modulation oscillator up to a pitch detector.  This is a quarter century ago, so I don't remember everything.

I s'pose you could derive that delay modulation signal from the input and likely LPF it before applying it to the instantaneous delay of the signal going out.


Detecting phase difference between a clean and a noisy signal

Here is how I would do it:

1) Find the frequency of the clean signal, this should be trivial in the time domain for a clean square wave.

2) Heavily smooth the noisy signal with a symmetric filter.  (IIR or FIR)

3) Calculate a single DFT bin for each signal using the length of a whole number of cycles for your frame.  This will make the harmonics leakless, and make your fundamental measurement more accurate.  The bin number you calculate should be the number of periods in your frame.

4) Calculate the phase for each signal relative to the frame at the fundamental frequency.  (Straight from the DFT bin value)

5) Take the difference in the two phases and translate if from an angle to a time using the frequency.

This is going to be way more efficient than some cross-correlation type calculation.


In Search of Voice Activity Detection Algorithms - What happens to the frequency or pitch or speech features as speaker moves away from microphone

The intrinsic fundamental frequency (pitch) of a speech signal do not change when the speaker moves away from the microphone (unless it moves away quite fast in which case it's a case of Doppler effect!) - it would be strange if we heard people's voice transposed up or down as they move away or as ambient noise increases!

However, as the distance from the microphone increases, the signal to noise ratio will increase, and this might affect the reliability of a naive pitch detector, so you might indeed get different result when estimating the f0 of a speech signal and the same signal with more noise. It doesn't mean that the f0 is different, just that the noise makes the estimate unreliable.

I am not aware of any naive feature which allows noise-robust speech/non-speech discrimination with a simple threshold rule. I think supervised machine learning would be the way to go for your problem.


In case of using a Doppler radar, do we have to use IQ demodulation to find the direction of a target?

Radar designer here: It sounds like you’re talking about  pulse-Doppler (PD) radar systems. For PD radars, the process is essentially as you described:

Generate a waveform (typically at IF) and then convert to RF
Transmit the waveform at RF
Receive the waveform at RF, and eventually mix it down to IF.
Apply IQ demodulation (for digital receivers a Hilbert transformer is typically used but that’s not an important detail). In this step, the received signal is added with a version that is 90 degrees out of phase with itself multiplied by the square root of -1, or the imaginary number. In general, you’ll have the generic signal s(t) = I(t) +jQ(t), where j is the symbol for the imaginary number. You can accomplish IQ demodulation in analog hardware, or you can use a digital receiver to accomplish it as well when you sample (via a Hilbert transform)
Collect the individual pulses and “stack” them in an array
Perform pulse compression and Doppler processing. Doppler processing is accomplished by computing the DFT across each pulse.

Since each waveform has phase which is recoverable with IQ data, the cross pulse DFT basically allows you to see the rate of change of that phase information, which is naturally the Doppler frequency. By being complex, yes you can tell whether a target has a negative or positive velocity, which is naturally very useful and pretty crucial for most tracking radars.

As for your question about needing to modulate the waveform: if by modulation you mean something other than a pure sinusoid, yes generally that’s a good idea for pulsed Doppler radars, as you’ll need good pulse compression. LFMs (linear frequency modulated) signals are a popular choice due to their nice properties, though other modulation schemes exist (NLFM, PSK, FSK, etc). If you just use a pure sinusoid, you effectively get no pulse compression, and you’ll have a pretty poor estimate of where targets may be, and the performance would be awful if you have several closely spaced targets.


There are several ways.

Start with taking the iFFT of $ S^{(1)} $ and $ S^{(2)} $.  Let's call the results $ m^{(1)} $ and $ m^{(2)} $, since they are modified.  They are back in the time domain and you wish to recombine them into one signal again.  Ideally, they would line up perfectly and you could just truncate one and keep the other.  This won't happen, so the trick is to fade the latter into the former over the overlap region.  At the beginning of the region the weighting of $ m^{(1)} $ should be one and the weighting of $ m^{(2)} $ should be zero.  At the end of the region, the reverse should be true.  The simplest solution is a linear transition, but that may not be the best.

For a 50% overlap, there is a nice property of a VonHann window in that the appropriate weighting is implicitely correct in the transition region and you can simply sum the two modified signals over the overlap region.  For other overlap sizes the same function can still be used, but it needs to be resized appropriately.  The window function can be applied in the time domain before you take the initial FFT, or after the iFFT is done.  Of course, which you do greatly affects the bin values, so it depends on the nature of your "touching" on which way is better.  The equivalent of a VonHann window can also be done in the frequency domain by calculating a new value for each bin by subtracting the average value of its neighbors and taking half that value.  The equation would be $Y[k] = -.25 * X[k-1] + .5 * X[k] - .25 * X[k+1]$ where the $X$s are the original bin values.  Doing this before or after your "touching" is the same as if you were windowing the signal before going into the frequency domain or after coming out.

The less touching you do, the more similar the two modified signals are going to be so the smoother the fading transition.  For too much touching, the mismatch between the two modified signals may be too great and the results won't be so good.  However, with a fading over the overlap region, the results are always going to be a smooth signal.

The VonHann window is good for spectrogram displays because it spreads the value of a sharp peak (near integer frequencies) to the neighboring bins and reduces the leakage of in between frequencies, so when you do a tone sweep through the frequency range, the display remains somewhat consistent and less frame size dependent.  Chances are that this same property will make it desirable for you to do the window function before you do your touching.

Hope this helps,

Ced

=================================================

Followup:

The process works like this:

1) Divide your signal into overlapping chunks

2) For each chunk

A) Apply Window function

B) Do the FFT

C) Apply touches

D) Do the iFFT

E) Reverse Window function

3) Recombine chunks into a new signal

Generally speaking, when you recombine the signal you will want a weighting function that goes from 1 to 0 across your overlap region so you are fading out of one chunk as you are fading into the next chuck.

If you use an overlap of greater than 50%, you will have regions where more than two chunks exist for the same time domain and you'll have to figure out an appropriate recombining strategy.  My understanding is that overlaps over 50% are in applications where the end result is the frequency domain and your interest is tracking frequency components.

You can use any window function you want, and any overlap amount.  If you follow the steps above, if you don't do any "touching" in the frequency domain, you should get your original signal back.  The advantage of the VonHann window using a 50% overlap is that the window function actually works as a weighting function so steps 2E and 3 are significantly simplified.

Check out this article:  "Overlap_Add_OLA_STFT_Processing"

The Hamming window is covered below figure 8.10.


Finding edges in a color image can be done by decomposing the image into its channels, finding the gradients separately and fusing them somehow. However, such approach doesn't incorporate the color components in a joint model. Luckily, there is a better way to do this, which is the structure tensor representation.

The color structure tensor describes the bi-dimensional ﬁrst order differential structure at a certain point in the image. It is specified by:

$
   S=
  \left[ {\begin{array}{cc}
   R_x^2+G_x^2+B_x^2 & R_xR_y+G_xG_y+B_xB_y \\
   R_xR_y+G_xG_y+B_xB_y & R_y^2+G_y^2+B_y^2 \\
  \end{array} } \right]
$

where subscripts denote the partial derivatives. This is a more precise description of the local gradients. Eigen-decomposition is then applied to the structure tensor matrix $S$ to form the eigenvalues and eigenvectors. The larger eigenvalue shows the strength of the local image edges (gradient magnitude) and the corresponding eigenvector points across the edge (in gradient direction). In other words, eigenvalues encode the gradient magnitude while the eigenvectors contain the gradient orientation information. Note that the eigenvectors and eigenvalues can be computed analytically for a structure tensor.

There are many other usecases of the structure tensor, so it is good to know. Last but not least, here is a MATLAB code, which uses this information to compute the gradients:

http://www.mathworks.com/matlabcentral/fileexchange/28114-fast-edges-of-a-color-image-actual-color-not-converting-to-grayscale/content/coloredges.m

After obtaining the gradients, you could further carry on your edge detection in a standard fashion (canny etc.).


Eigenfunctions of discrete-time LTI systems

i think the only family of functions that are eigenfunctions to LTI systems the sole exponential.

for continuous-time, if it's LTI:

$$ y(t) = \int\limits_{-\infty}^{+\infty} h(u) \ x(t-u) \ du $$

the eigenfunction is $x(t) = A e^{st}$ and the output is

$$ y(t) = H(s) \ x(t) $$

where

$$ H(s) = \int\limits_{-\infty}^{+\infty} h(t) \ e^{-st} \ dt $$

for discrete-time LTI systems, it is similar

$$ y[n] = \sum\limits_{m=-\infty}^{+\infty} h[m] \ x[n-m] $$

the eigenfunction is $x[n] = A z^n$ and the output is

$$ y[n] = H(z) \ x[n] $$

where

$$ H(z) = \sum\limits_{n=-\infty}^{+\infty} h[n] \ z^{-n} $$

often $s=j\Omega$ or $z=e^{j\omega}$ but it wouldn't have to be that.


Measuring similarity between two very similar pictures

This measure really depends on how your images look like. A very basic thing would be the calculate the point-wise difference between both images and summing them up. The smaller this value is, the more equal they are. Possibly, you can perform some Gaussian Smoothing and contrast equalization before. But, that's really depending on the images.


Consider a continuous-time periodic signal $x_1(t)$ whose fundamental period is $T_1$, fundamental radian frequency is $\omega_1 = \frac{2\pi}{T_1}$ and CTFS (continuous-time Fourier series) coefficients are: $$ x_1 \longleftrightarrow c_k = \frac{1}{T_1} \int_{<T_1>} x_1(t) e^{-j\frac{2\pi k}{T_1}t} dt ~~ , ~~\text{for}~~ k=0,\pm 1, \pm 2...$$

We want to find the new CTFS coefficients, denoted as $d_k$, associated with the new signal $x_2(t)$ , with period $T_2$ and fundamental radian frequency $\omega_2 = \frac{2\pi}{T_2}$, and related to $x_1(t)$ as follows:
$$x_2(t) = x_1( a(t-b) )$$ where $a$ and $b$ are constants, assuming $a > 0$. Then we say ; $$ x_2 \longleftrightarrow d_k $$

We shall proceed from the above definition to show the relation between the coefficients $d_k$ and $c_k$ as follows:
$$ d_k = \frac{1}{T_2} \int_{-T_2/2}^{T_2/2} x_2(t) e^{-j\frac{2\pi}{T_2} k t} dt $$

First note that $T_2 = T_1/a$. And replacing $x_2$ with $x_1$ yields the following:

$$ d_k = \frac{a}{T_1} \int_{t=-T_1/{2a}}^{t=T_1/{2a}} x_1(a(t-b)) e^{-j\frac{2\pi a}{T_1} k t} dt $$
Now apply a substitution $t-b = t'$, and then replace $t'$ with $t$ again (note that $dt'=dt$)

$$ d_k = \frac{a}{T_1} \int_{t=-\frac{T_1}{2a}-b}^{t=\frac{T_1}{2a}-b} x_1(at) e^{-j\frac{2\pi a}{T_1} k (t+b)} dt $$ We know that for a periodic signal, as long as the integration involves a single period, we can redefine its limits:

$$ d_k = e^{-j \frac{2\pi a}{T_1} k b} \left( \frac{a}{T_1} \int_{t=-T_1/{2a}}^{t=T_1/{2a}} x_1(at) e^{-j\frac{2\pi a}{T_1} k t} dt \right)$$
Now apply substitution $at = t'$ (note that $dt' = a dt$) and then replace $t'$ with $t$ again which yields:

$$ d_k = e^{-j \frac{2\pi a}{T_1} k b} \left( \frac{1}{T_1} \int_{t=-T_1/{2}}^{t=T_1/{2}} x_1(t) e^{-j\frac{2\pi}{T_1} k t} dt \right) $$
Finally recognise that the integral equals $T_1 c_k$ and the relationship between $d_k$ and $c_k$ is:

$$ d_k = e^{-j \frac{2\pi}{T_1} a b k} c_k $$

Note that the coefficients $d_k$ are for a periodic signal with fundamental frequency $\omega_2 = \frac{2\pi}{T_2}$, whereas $c_k$ are for the frequency $\omega_1 = \frac{2\pi}{T_1}$. So this means that eventhough, for example, $|d_3|=|c_3|$, i.e., both coefficients have the same magnitude for $k=3$, those coefficients are placed at different frequencies; $d_3$ is associated with $\omega = \frac{2\pi}{T_1}6$ whereas $c_3$ is associated with the frequency $\omega=\frac{2\pi}{T_1}3$... It can be seen that the new coefficients $d_k$ have the same magnitudes, i.e., $|d_k| = |c_k|$, but different phases.


STFT output frequency does not match sin wave frequency

You have to pass sr=fs as well as hop_length=hop_length to the specshow function, as its default values are different from yours.


The algorithm is basically based on the assumption that noise is "white", ie. in the expectation value, the noise energy in every DFT bin is equal. Under that assumption, and assuming a "long enough" observation, there has to be the same noise floor in every bin.

Sound on the other hand, is a highly correlated signal. So the assumption is that there's at least one DFT bin at a time where there's no significant sound signal.

As you might imagine, that's a bit of an oversimplification on both sides: There's a lot of music that actually uses Noise, or broadband signals, and audio noise typically is everything but white.

In essence: this boils down to the good ol' "What is my SNR?" question:

You, as the DSP designer, will have to model your signal, and your noise. After that, the question kind of answers itself. It's really more of a question of semantics than engineering: For example, take psychedelic rock/stoner/doom metal: wideband noise in the string instrument's track is very essential to that genre; is it noise, then? And if you get a recording of a classical guitar piece, wouldn't you consider the same noise in the signal part of the noise or music energy?


The raw data you get from the camera is nearly linear in the colorspace spanned by the camera CFA filters. The 8-bit/sample monitor you have assumes the image is gamma-compressed, most likely within the sRGB colorspace. These two colorspaces are different, therefore you need to map between the two in some way. Understanding this, lets address each of your issues:

  The RGB image looks is almost completely black.

The maximum 10-bit value is 1023, the maximum 16-bit value is 65535. To Convert 10-bit linear to 16-bit linear you multiply by 65535./1023 and round to nearest integer. Alternatively it might be easier to work with floating point directly, so just divide by 1023 to map your data to the [0,1] range.

  the RGB image looks ... rather green

It is common for the camera sensors to be more sensitive to green, so the green response is higher than red and blue, giving an overall green tint. To get physically correct measurements of the incoming light within a known colorspace (say the ITU-R BT.709 primaries used by your monitor) you need to figure out the color transformation matrix from the camera colorspace to your preferred colorspace. There are different ways to accomplish this: get a color-calibration chart, contact the sensor manufacturer, or measure the exact spectral response using other means (like a monochromator).

That being said, if you don't seek for physically correct measurements but rather just for a good-enough picture, you might skip the above calibration step and instead do only a white-balancing step: find an ought-to-be gray object in the frame and divide each pixel by the RGB of that object, thus cancelling the green tint. This on its own, without a color transform, will not reproduce colors exactly though.

  The overall intensity is very low

This is due to the linear versus gamma discrepancy. You do all your math in linear color space with high bitness, and as a last step, convert it to 8-bit sRGB for display. See sRGB for the exact formula.


How can I automatically find the boundary of a high intensity region inside an image?

One method that I have had lots of success with, and others have too, is Circular Shortest Path (CSP). Here is one paper, but there are quite a few more and possibly some code floating around the web. Another technique that are perhaps more popular, and for which code might be more common, are snakes (active contours).

These techniques work by finding a connected path that minimises a certain function. In your case, it would be to find the path that goes through many edge pixels while still remaining somewhat round. The roundness constraint is what keeps the path from going off down the narrow gaps. Snakes are different from CSP in that the boundary is 'grown' out from a seed area; CSP performs a polar transformation about the centre point.

Another possibility is some kind of morphological operator. These can be used to find regions of a particular size. One would convert the image to binary, the apply a morphological opening using a structuring element of a big enough size to fit inside the area of interest, but not anywhere else.


No, you're right. It wants to know how many samples you have for each symbol in the BPSK stream. This is the ratio of the sample rate to the symbol rate.


Why zero padding the 2-d DFT interpolates images in spatial domain?

For the sake of simplicity, I'll explain the 1-D case; the 2-D case is completely analogous. Let $x[n]$ be a finite length sequence with $n\in[0,N-1]$. Its discrete Fourier transform (DFT) is

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi}{N} kn}\tag{1}$$

The sequence $x[n]$ can be obtained from $X[k]$ via the inverse DFT (IDFT):

$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j\frac{2\pi}{N}kn}\tag{2}$$

Note that the sequences defined by $(1)$ and $(2)$ are periodic with period $N$. Let's use $\tilde{X}[k]$ to denote the periodic sequence of DFT coefficients. For odd $N$, $\tilde{X}[k]$ just equals (the periodic continuation of) $X[k]$. For even $N$, we split the coefficient $X[N/2]$ in two equal terms $\frac12 X[N/2]$, such that for even $N$, the periodic sequence $\tilde{X}[k]$ is defined as

$$\tilde{X}[k]=\\\big[\ldots ,\underbrace{X[0],X[1],\ldots,\frac12 X[N/2],\frac12 X[N/2],\ldots,X[N-1]}_{N+1\textrm{ elements}},\ldots\big]\tag{3}$$

With that definition, the period of $\tilde{X}[k]$ is always odd. Let's denote that odd period by $2K+1$. For odd $N$ we have $K=(N-1)/2$, and for even $N$ we get $K=N/2$.

Using $\tilde{X}[k]$ and $K$ as defined above, we can rewrite $(2)$ as

$$x[n]=\frac{1}{N}\sum_{k=-K}^{K}\tilde{X}[k]e^{j\frac{2\pi}{N}kn}\tag{4}$$

With these preparations we can now interpret the coefficients $\tilde{X}[k]$ as Fourier coefficients of a continuous periodic function $x_c(t)$ with period $NT_s$, where $T_s$ is an arbitrary sampling interval:

$$x_c(t)=\frac{1}{N}\sum_{k=-K}^{K}\tilde{X}[k]e^{j\frac{2\pi}{NT_s}kt}\tag{5}$$

The sequence $x[n]$ can be obtained by sampling $x_c(t)$:

$$x[n]=x_c(nT_s)\tag{6}$$

But we can also obtain other, more densely sampled sequences by sampling $x_c(t)$. We choose the new sampling interval such that the new sequence is also periodic with period $M>N$:

$$\hat{x}[m]=x_c\left(m\frac{NT_s}{M}\right)=\frac{1}{N}\sum_{k=-K}^{K}\tilde{X}[k]e^{j\frac{2\pi}{M}km}\tag{7}$$

The new length $M$ sequence $\hat{x}[m]$ can be interpreted as an interpolated version of the original length $N$ sequence $x[n]$.

Note that $(7)$ can be written as a length $M$ IDFT with a zero-padded sequence $\tilde{X}_0[k]$:

$$\begin{align}\hat{x}[m]&=\frac{M}{N}\left(\frac{1}{M}\sum_{k=0}^{M-1}\tilde{X}_0[k]e^{j\frac{2\pi}{M}km}\right)\\&=\frac{M}{N}\cdot \textrm{IDFT}_M\big\{\tilde{X}_0[k]\big\}[m]\tag{8}\end{align}$$

with

$$\tilde{X}_0[k]=\big[\tilde{X}[0],\ldots,\tilde{X}[K],\underbrace{0,\dots,0}_{M-2K-1\textrm{ zeros}},\tilde{X}[K+1],\ldots,\tilde{X}[2K]\big]\tag{9}$$

In terms of the original DFT coefficients $X[k]$, the zero-padded sequence $\tilde{X}_0[k]$ is given by

$$\tilde{X}_0[k]=\big[X[0],\ldots,X[(N-1)/2],\underbrace{0,\dots,0}_{M-N\textrm{ zeros}},\\\qquad X[(N+1)/2],\ldots,X[N-1]\big]\tag{10}$$

for odd $N$, and by

$$\tilde{X}_0[k]=\big[X[0],\ldots,\frac12 X[N/2],\underbrace{0,\dots,0}_{M-N-1\textrm{ zeros}},\\\qquad \frac12 X[N/2],\ldots,X[N-1]\big]\tag{11}$$

for even $N$.


Image Histogram - What Is the Interpretation of PDF

Well, If you model your image as a realization of a random variable generator then the Histogram is the best estimation (Assuming no other information like prior, etc..) you have for the PDF of the random variable.

For instance, you can see this model is used when doing Histogram Equalization (Transforming the realization into a realization of Uniformity Generator).

Pay attention that this is a very simple model.
For instance it doesn't take care of the correlation between adjacent pixels in image.

Indeed your interpretation is correct given the model.


Say you have two pulses, $p(t)$ and $q(t)$, such that the following three statments are true:

\begin{align}
\int_{-\infty}^\infty p(t) p(t) dt &= 1,\\
\int_{-\infty}^\infty q(t) q(t) dt &= 1,\\
\int_{-\infty}^\infty p(t) q(t) dt &= 0.
\end{align}

Then these pulses are said to be orthonormal. Now, let us say you want to transmit two information-bearing symbols $a_1$ and $a_2$. (The symbols carry information when the receiver does not know their values in advance). You can acomplish this by transmitting the signal $s(t) = a_1 p(t) + a_2 q(t)$.

How does the receiver recover the symbols? By performing these operations:

\begin{align}
a_1 &= \int_{-\infty}^\infty s(t) p(t) dt, \\
a_2 &= \int_{-\infty}^\infty s(t) q(t) dt. \\
\end{align}

I think it's clear that this scheme can be extended to any number of symbols; all you need to do is find a set of orthonormal pulses of the required size. The most common approach is to find a "prototype pulse" $p(t)$ with the property that, when time-shifted $lT_p$ seconds with $l$ integer, it is orthonormal to all other pulses time-shifted $mT_p$ seconds, for any integer $l \neq m$.  So, you can transmit any number of symbols with the signal

\begin{align}
s(t) &= a_1 p(t - T_p) + a_2 p(t - 2T_p) + a_3 p(t - 3T_p) + \cdots \\
&= \sum_k a_k p(t-kT_p).
\end{align}

The symbols $a_k$ are recovered in the same way as before, but this operation can be implemented with a filter "matched" to the prototype pulse $p(t)$ (see numerous questions related to matched filters on this site). When the matched filter's output is sampled at time $t=kT_p$, you obtain back the symbol $a_k$.

Note that the filter's output is not necessarily sampled at a rate that meets Nyquist's sampling theorem (I have said nothing about the bandwidth of $p(t)$ or $s(t)$), but that is OK because you're not interested in reconstructing $s(t)$ from its samples, all you want is the symbols $a_k$.


Given input and output time series, obtaining filter design

If your system is linear, then you can find the impulse response $h(n)$ of the system using the convolution sum $y(n)=\sum_{k=0}^{N} h(k)x(n-k)$, which can be written in matrix-vector form as $y=Xh$, where $X$ is a convolution matrix of the input data and $N$ is the order of the system, say $N=10$ for example. Since you know $X$ and $y$, you can find $h$ via the Pseudo inverse, assuming you have created an over-determined system $y=Ah$. The pseudo-inverse is $h=(X^T X)^{-1}X^T y$. After finding $h$, find the FFT of $h$ or take the Laplace transform of $h$ to get the transfer function. I hope this helps you.