Well, for a signal $x[k]$, the energy is:
\begin{equation}
E=dt*\sum_{k=-\infty}^{\infty}x^2[k]
\end{equation}

You want the energy to be one by dividing $x[k]$ through a value $a$, hence:
\begin{equation}
1=dt*\sum_{k=-\infty}^{\infty}(\frac{x}{a})^2[k]
\end{equation}
which by multiplying with $a^2$ becomes
\begin{equation}
a^2=dt*\sum_{k=-\infty}^{\infty}x^2[k]
\end{equation}
Comparing this to the first equation, it can be easily seen that $a^2=E$ and since $E$ is always positive, we can conclude $a=\sqrt{E}$. And that is why you divide by the square root of $E$.

From there it should be easy to see why multiplying this with $\sqrt{E'}$ makes a signal of energy $E'$.


Gibbs phenomenon in Hamming's digital filters

I think the main issue is you are jumping ahead of yourself.  You probably remember or read somewhere that the Fourier transform of a $rect$ function is a $sinc$ function.  This is true; however, no where in this section does he mention Fourier transform!  In fact, what he is doing is not Fourier transform.

What he does in this section is to represent any periodic function as a Fourier series:

(1)    $f(\theta)=\frac{a_0}{2}+\sum_{k=0}^\infty{(a_k \cos k\theta + b_k \sin k\theta)}$

The key here is that this function doesn't contain every frequency.  It contains frequency 0 (the DC) and $\frac{k}{2\pi}$ where $k=1,2,3,...\infty$.  This is actually very important.  This function is always periodic with a period of $2\pi$.  The coefficients $a_k$ and $b_k$ are almost sampled values of the Fourier transform (but not quite, why? homework exercise! :p).  You have countably many of these coefficients.  Later in the book, you will learn that when you sample data in one space (frequency or time) it necessarily makes the counter part in the other space (time or frequency) periodic.

In the first case, he does the Lanczos smoothing derivation where he averages the function by running a rectangular window through it (convolving with rect).  What he shows is, not surprisingly, that the coefficients get multiplied by this term:

(2)    $\frac{\sin k\pi/N}{k\pi N}$

which should look very familiar to you, of course, because it is the $sinc$ function.  However, what you are missing is that $k$ is discrete!  It is actually a sampled version of the $sinc$ function.

Effectively, he convolves a function with a $rect$ and shows that the coefficients of the resulting Fourier series (read loosely as: sampled Fourier transform) is a sampled $sinc$ function.  No surprise there.  Convolution thm says convolution in time turns into multiplication in frequency.  Fourier transform, which you know, of $rect$ is $sinc$, so convolution by $rect$ is multiplication by $sinc$ in frequency space.

In the next section, he does something different.  He takes the Fourier series (read: sampled Fourier transform) and removes all the higher frequency coefficients.  In effect, he is taking the Fourier transform, multiplies by a $rect$, and then samples it.  For simplicity, he sets all the Fourier coefficients that did not get discarded to $1$.

What he's left with is this:

(3)    $h(\theta)=\frac{sin(N+1/2)\theta}{sin(\theta/2)}$

And you ask, why isn't this a $sinc$ function?  Can you answer it now?

The quick answer is because what is applied in frequency domain is not just a truncation, it's a truncation and a sampling operator.  What you know is that when you truncate (i.e. multiply by rect) in frequency domain, the time domain gets convolved with a $sinc$ (by Convolution thm and Fourier transform of $rect$), but this is without sampling.

As for why the formula looks the way it does, there are two ways to look at it.  The first, which he shows, is that you can just sum the Fourier series from $-N$ to $N$, and that's what you get.

The second, which is more profound and may come up later, is that (*) when you sample a function in one space (say frequency), the corresponding function in the other space (say time) becomes the sum of shifted versions of itself.  In fact, it probably won't come up exactly like that.  The typical scenario is that sampling time domain creates periodic replication in the frequency domain (btw: this is the reason for what people call aliasing).  However, you can apply the duality property of Fourier transform to get (*).

Does equation (3) make sense now?  It is periodic.  The closer you get to 0, the more it looks like a $sinc$ function.

So, an exercise for you is to derive equation (3) by sampling and truncating in frequency space and applying inverse transform.


How to Extract High Frequency and Low Frequency Component Using Bilateral Filter?

Similar to one dimensional signals, low frequencies in images mean pixel values that are changing slowly over space, while high frequency content means pixel values that are rapidly changing in space.

For example, the following image has strong low frequency components: You can intuitively see how I simply have a sin-wave propagating at some low frequency.

In contrast, this image below is composed of a sin wave at double the above frequency. Notice how this manifests itself on the image:

Any image can have any number of low and high frequency components together as well. For example, and image like this has both low and high frequency components:

You can see how you have a low-frequency 'trend', but also a lot of high frequency detail across the image. (Roughly speaking, in images, sharp transients like edges correspond to high frequencies, while long unchanging spaces correspond to low frequencies).

So now, the point of an application of a bilateral filter (which is simply a convolution of your image with a gaussian kernel), is to remove high frequency components, and retain your lower frequency components. So in this case, what happens if we convolve the above image, with a gaussian (bilateral) filter that looks like this?

If we use this as our kernel, we will get the following result:

In other words, we have removed the high frequency components of the original image, but kept the low frequency components, since the variance of the gaussian kernel was chosen appropriately.


Is ergodic in mean a property defined only for WSS stochastic processes?

Time average is defined over a time interval, T which will always result in a constant value. Keeping this constraint in mind, ensemble mean is always forced to be constant for satisfying the condition of ergodicity.


Find Preamble (pattern) in signal, using Corelation

You should be able to see an peak correlation as from your plot and the length of your pattern it does not appear to be noisy enough such that you would not see a peak (unless your noise is completely in the time domain as in excessive jitter that I cannot visibly see).

I suspect that the pattern you are using is extending far beyond the length of the preamble such that your data payload is corrupting the final cross correlation result (or you are not using just the preamble pattern to correlate with your data). It goes without saying that the preamble and data need to be running at a similar rate (although you can have a significant clock offset and still correlate; correlation amplitude versus frequency offset is a Sinc function with the first null at a clock offset of 1/T where T the length of your correlation sequence; the longer the sequence the more accurate your clocks need to be). From your result if I read the numbers on the plots correctly your correlation result is close to 60 on average, and from your description you are correlating over 32 symbols (4x 0xAA).  If the sequences were not correlated, and the average for your data and pattern are 1/2, then this suggests your sampling rate is approximately 8 samples per symbol if you are indeed correlating over the full 32 symbols.

You will also benefit by removing the mean from your data (such that it is bipolar in the correlation computation) as you will see in my example below. Also as shown this does not appear to be the reason you are not getting any result (as shown by my result even with the mean present).

Trying to duplicate your conditions, here is the result I got when I correlated the 4x 0xAA pattern to a random sequence of 0's and 1's for 5000 samples: (although I reversed the correlation from what you did such that the preamble correlation would appear first):

And then with the preamble at the beginning and repeated in the middle of the random sequence (note the noise is very big in this case; as I added the random data to the preamble, just to show what you should expect to see in extremely noisy conditions:

And by removing the mean of the signals prior to correlation the detection becomes quite easy even in the presence of the very significant noise (the noise is a random pattern of 0's and 1's at sample rate):

Also for completion, below shows the result in the case of correlating over just a portion of your sequence (1x 0xAA) with the same preamble in the sequence (4x 0xAA); 20Log(4) less processing gain but showing that a clear detection is still possible with the amount of noise I had added (and such an approach may be done if the rest of the preamble is reserved for timing recovery; meaning by design we need to detect early in the preamble):

Not directly related to why you are having issue, but a challenge also to be aware of with this pattern is that while great for timing recovery (as there is a transition at your data rate on every symbol), your preamble does not have a unique autocorrelation as you can see by the cyclical pattern in the plots above.

The same processing gain in SNR is achieved as in other patterns at the position in cross-correlation where the delay between the patterns is zero , $R_{xy}(0)$, (the gain in SNR is 20Log(N) where N is the number of samples where the noise is still uncorrelated; usually due to symbol rate filtering this is the number of symbols) and if you use a circular correlation you can also achieve that SNR advantage at every offset where the sequence aligns as described above as well. For timing critical applications where it is important to know the actual position within the sequence, Barker codes or other similar pseudo-random sequences are used that have the property that correlation is maximum at shift = 0 and close to 0 everywhere else (such as maximum length sequences).  An example 11 bit Barker code is 10110111000.  If you treat the 0's as -1 you will see that this sequence will add to +11 when they are added and in a circular correlation will add to +1 everywhere else. Note however the challenge for clock recovery in such a sequence given the longer runs of ones and zeros compared to 1010101010.


You are apparently in the context of no-reference, reference-free or blind image quality assessment. The topic is quite active, and I am not sure people have already a completely accepted framework for that. Multiple distortions may affect images: random noise, compression artifact, static blur, motion blur, etc. They require different metrics (benchmark data here). Some references to start with, and help you for more focused questions:

No-Reference Image Quality Assessment in the Spatial Domain, 2012
Making a “Completely Blind” Image Quality Analyzer, 2013
No-reference image quality assessment based on spatial and spectral entropies, 2014
MDID: A multiply distorted image database for image quality assessment, 2017

The first three papers come from the same group, an hint that the paint is still wet.


Clarification Needed: Power Spectral Density of Real Signal

This stems from Parseval's theorem. For any signal, $\int_{-\infty}^{+\infty}|x(t)|^{2}\mathrm dt = \int_{-\infty}^{+\infty}|X(f)|^{2}\mathrm df$ and hence for a real signal, the one sided PSD must be multiplied by a factor of 2 because of Symmetry. Intuitively, the Fourier Analysis of any signal provides us with a spectrum which has Negative frequencies too because of the fact that Sines and Cosines can be constructed using Euler's theorem of using Complex Exponentials of positive and negative frequencies; and it so happens for any real signal that all you need are complex exponentials which are conjugates of each other to construct them.

The negative portion is not exactly "redundant". The signal still owns those "negative" frequencies and hence anything present there accounts to its energy.


When you slide a rectangular signal over another rectangular signal and take the area under the curve, i.e., perform convolution, at first there is no overlap. Then, the overlap gradually starts increasing in a linear fashion. When the two rectangular signals are exactly on top of each other, that is the point of maximum overlap. After than, there is a gradual decrease until there is no overlap. The resulting convolution signal is then a triangular signal.

Update: A rectangular signal is given by
\begin{equation*}
\prod (t) =
            \begin{cases}
1 &  0 \le t \le 1 \\
0 & otherwise
            \end{cases}
    \end{equation*}
Its convolution is given by
$$z(t) = \int _{\tau=-\infty}^{\infty} \prod(\tau) * \prod(t-\tau) d\tau$$
which turns out to be
\begin{equation}
z(t) = \begin{cases}
0 & -\infty \le t \le 0 \\
t & 0 \le t \le 1 \\
2 - t & 1 \le t \le 2 \\
0 & 2 \le t \le \infty
\end{cases}
= \Lambda(t)
\end{equation}


Gaussian noise generator algorithm trouble

I agree with you, this looks like a flawed description of an inverse transform generator.

The $p[i]$ shouldn't be like described, but should be the prefix sum of these, corresponding to the cumulative distribution.
The $\arg\min$ should indeed be taken on the absolute value (or allow positive values only).
As the pdf is computed on the positive $i$ only, the range of the uniform drawing should be $[0,\frac12]$ (or the $p[i]$ scaled by $2$.)

(When the $p[i]$ are not cumulated, they are in $[0,\frac1{\sigma\sqrt2\pi}]$, and the range $[0,1]$ doesn't make sense.)

For better efficiency, you can precompute a lookup-table of the $i$ indexes for sufficiently dense $q_1$.


It's not correct to say that both expressions become equal when the signal is periodic. What is correct is that - obviously - everything one can know about an $N$-periodic signal is represented by the expression on the left of your equation, and consequently, both expressions are equivalent in the sense that one can be obtained from the other.

Note that that the infinite sum on the right-hand side does not converge in the conventional sense if the signal is periodic. The Discrete-Time Fourier Transform (DTFT) - that is what the right-hand side represents - of an $N$-periodic signal can be expressed by using Dirac delta impulses:

$$X(\omega)=\frac{2\pi}{N}\sum_{k=-\infty}^{\infty}X[k]\delta\left(\omega-\frac{2\pi k}{N}\right)\tag{1}$$

The coefficients $X[k]$ can be obtained from the expression on the left-hand side of your equation be setting $\omega_k=2\pi k/N$. In that sense, both expressions convey the same information.


How to efficiently calculate points of intersection of a straight line and a contour?

Calculate the straight lines angles with respect to the x-axis

Calculate the angles of each and every point on the contour

Search for points on the contour with same angles as the lines with respect to the origin.

These are the points of line-contour crossing.


How to use clustering file after training

You can use cv::FileStorage. You can follow this post here and this documentation here.


How to get Step Response from an Impulse Response in Simulink?

Your system H(s) is continuous

$$ H(s) = \frac{s}{s+1} $$

Using a discrete impulse does not make a lot of sense. I suspect Simulink added a "zero-order hold" or ZOH block between your discrete impulse and your system $H(s)$. As such, your impulse response will be affected by the sampling period of your discrete impulse source. If you decrease the sampling period, the second curve shoud match the step response more closely, although they will never be exactly the same.

Edit :

Decrease the width of your impulse to say 0.1 second and increase the amplitude to 10 so that the area under the curve stays 1.


Design of a low pass filter with Red Pitaya's FPGA

There are 2 steps for the solution
Step 1 - Design some kind of low-pass filter to measure the DC component.
In order to perform step 1 you need to answer a few questions :
What do you mean by DC? For example, measure the DC for the last 10 seconds?
What is the input sampling rate? What is the desired output sampling rate?
You need to provide this information, we cannot guess it for you.
Step 2 - Implement the filter in Verilog.
There are examples on opencores. I can also give you a hand if you need help with implementing a FIR or IIR filter in Verilog.
You could also look at the Xilinx FIR Cores in the Core generator.


Why do third-octave bands have a ratio of cube root of two?

A fixed interval corresponds to a fixed ratio of frequencies (at least in equal temperament, which is the most common tuning system in Western music). One octave is a factor of $2$. One semi-tone is a factor of $\sqrt[12]{2}$ because there are $12$ semi-tones in one octave. So one third of an octave (i.e., $4$ semi-tones) corresponds to a factor of $\sqrt[3]{2}$ because $\sqrt[3]{2}\cdot \sqrt[3]{2}\cdot \sqrt[3]{2}=2$ (i.e., one octave).


Can you explain about the following formulas? in Signals and System

Your question is about the Fourier series decomposition of a function (called signal here). First, this does not apply to "any signal" in theory. In practice, it applies to most pratical "periodic" signals. Periodicity is quite important, althought Fourier development may exist for not strictly periodic functions.

In loose terms, it tells you that if a signal is periodic, then inside each period the signal can be expressed as a sum (a linear superposition) of the most simple periodic functions, sines and cosines, with frequencies taken as linear factors of the base frequency.

The $a_0/2$ term comes from the average of your signal over one period $T$. This is due to the fact that sines and cosines do not see constant signals: they are orthogonal to them, since the integral of a sine or a cosine over one period exactly sums to $0$.

The other terms $a_n$ and $b_n$ are related to orthogonality as well, and i guess that you can, with great profit, read the following Wikipedia: Fourier series page in your perferred language.


Reversing the order of up and downsampling

The problem with downsampling is that it can be lossy -- since you're reducing the sampling rate, you can introduce aliasing. So, you can reverse the order whenever downsampling does not result in aliasing.

For example, say your discrete-time signal $x[n]$ contains energy in frequencies up to $f_N/3$, where $f_N$ is the Nyquist frequency. Then, downsampling by a factor of 2 will not affect the signal: the "new" Nyquist frequency is $f_N/2$, so there is no aliasing.

See also Which order to perform downsampling and filtering?.


Why exactly is the Guassian filter a low-pass filter?

Because the Fourier transform of Gaussian in the spatial/temporal domain is also a Gaussian in the frequency domain. Since a Gaussian drops off rapidly with frequency, it's a low pass filter.


Each HRIR is meant to represent the acoustic impulse response from an audio source at a specified angle (relative to the head) to the ear canal.  There is no such thing as a 3D signal, but if you want to simulate audio from a given azimuth and elevation angle, just convolve your audio source with the corresponding HRIR as follows:

left_channel = conv(audio_signal,hrir_left_az_elev)
right_channel = conv(audio_signal,hrir_right_az_elev)

In this case, I am using hrir_left_az_elev to signify the HRIR that corresponds to the left-ear impulse response to a sound event at azimuth angle az and elevation angle elev.  How you actually code this up will depend completely upon the format of your HRIR data.


The Box–Muller method uses two independent random numbers distributed uniformly on $(0,1)$, then uses such two uniformly distributed numbers to produce values with normal distribution.

Zero-mean indicates the average value of the detected signal at each time point is zero. In digital signal processing, zero often refers to $2^{n-1}$ where $n$ is the bit number of digitizer. For example, if you have a 16-bit signal ($65536$), then $32768$ represent the zero in double, and those numbers that less than $32768$ are negative.


Extracting filter coefficients from raised-cosine frequency response

This is a well-known filter: it's called raised cosine filter. The best strategy is to explicitly compute the inverse Fourier transform of the given frequency response. The corresponding result can be found in the wikipedia article linked to above, and it's also given in this article, which has some additional information on pulse-shaping filters.

Since the corresponding impulse response is of infinite length, you need to truncate it symmetrically. Since the given filter order is $24$, you have 25 filter taps (assuming a common FIR solution). This means one tap at $n=0$, and $12$ taps for $n>0$ and $n<0$, respectively. For truncating the ideal impulse response, you can use any type of window. The most straightforward way is simple truncation, i.e., using a rectangular window.

You can use Matlab to evaluate the formula for the filter's impulse response, and for applying a window.


Convolution theorem for "weighted" convolution?

1) The terms autoconvolution or selfconvolution can be found in the literature, without the kernel $s$. So you would probably have a weighed, kernel or windowed self-convolution.
However, i foresee two potential caveats on your formulation. First, i do not understand why you require a two-dimensional periodic signal (and nothing on the kernel $s$?). This could result in infinite values of your integrals. Second, don't you assume some kind of symmetry in $s(r,r′)$?

2) Suppose you can solve your system with $s=1$, then with appropriate symmetry on $s$ you may be able to factorize $s$ and slit it over $g(r')$ and $g(r-r')$. Unfortunately the convolution theorem is likely to give you a family of solutions $g$ knowing $f$ only.

3) You may find useful to explore a multiscale form of your weighted self-correlation operator.

Affine invariant pattern recognition using Multiscale Autoconvolution
Rahtu E., Salo M., Heikkilä J.
IEEE Trans Pattern Anal Mach Intell. 2005 Jun;27(6):908-18.

An iterative deautoconvolution algorithm for nonnegative functions
Kerkil Choi and Aaron D Lanterman
2005 Inverse Problems 21 981

Autoconvolution and panorama: Augmenting second-order signal analysis
Douglas, S.C. and  Mandic, D.P., 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)


How to perform a phase recovery for higher order modulation schemes

This post provides many details on implementing carrier recovery for a higher order QAM system (or even QPSK).
High modulation index PSK - carrier recovery
Using a PI loop filter where in addition to a simple accumulator, a direct (proportional) path is also summed in results in a 2nd order Type 2 tracking loop as required to track phase offsets to zero error under static frequency offset conditions. Symbol-aided synchronization (where a known pilot is used) is recommended for higher order QAM to minimize acquisition time. The coefficients of the loop filter are set to stabilize the carrier recovery loop and also set the appropriate loop BW.
See this post for further considerations on setting the loop BW properly:
Loop bandwidth for symbol timing recovery


Why do I keep seeing this equation for $E_b/N_0$

It's just a matter of finding out the assumptions behind the equation. Sometimes, unfortunately, those assumptions are not made explicit.
For the first equation you present, the one you have issues with, the assumptions are (off the top of my head):

Orthogonal, ideal sinc pulses are used as pulse shape.
The signals involved are strictly baseband, or, if passband, then the upconversion is suppressed-carrier, and the entire process is ideal.
The transmitted signal is narrowband; there are no multicarriers.
The channel is AWGN; there is no signal distortion, not even attenuation.
There is no coding and the transmitted bits are uncorrelated.
The received signal is matched-filtered. The matched filter is ideal (i.e. an infinite sinc pulse).

Under these assumptions, that equation is correct.
If you change those assumptions, then the relationship changes, as you found. I didn't check every one of you equations carefully but in general it seems like you are on the right track.


How to apply the gain factor G to SOS coefficients when designing a IIR filter on sptool (MATLAB)?

The easiest way would be
y = prod(G)*sosfilt(SOS,x);

An alternative would be put to the cumulative gain into the first section.
sosScaled = SOS;
sosScaled(1,1:3) = sosScaled(1,1:3)*prod(G);
y = sosfilt(sosScaled,x);


The information about the frequency of the sine seems to be missing. This could be an hint for you to compute the DFT exactly (not with Matlab simulation), for a given frequency $f$ and a number $s$ of seconds. Since the sampling frequency $f_0$ gives a integer multiple of samples inside one second, you should be able to find manageable formulas.

As commented by @Juancho, there is a lack of homogeneity in your time base: a frequency times a duration in seconds is not homogeneous to a time (it is unitless). A flexible solution is:

fs=10; % sampling freq
lFreq = [3 3.5 4]; % list of frequencies
lDuration = [1 2 4 10]'; % list of durations

for iDuration = 1:length(lDuration)
    timeBase = [0:1/fs:lDuration(iDuration)]';
    dataCosine = cos(timeBase*lFreq*2*pi); % Three cosines
    data = sum(dataCosine,2); % sum the three cosine components
    figure(iDuration)
    plot(timeBase,data,'x-');
    ylabel('Amplitude (a.u.)')
    xlabel('Time (seconds)')
    title(['Duration: ',num2str(lDuration(iDuration)),' s'])
    grid on
end


Time frequency analysis for a spreading signal

The overlap is not important, it is just necessary that the spreading sequences  are pseudo-orthogonal.

Can someone confirm?


Finding Signals when the baseline varies greatly by signal set

This sounds like a job for a simple high-pass filter, followed by a simple threshold, to me.

1/3 billion records isn't much data if you're used to radio signal sampling rates, so computationally, this should be most definitely doable.

So, with scipy.signal there's a reasonable signal processing library for Python, but I'm from the GNU Radio world, where you can use Python to plug together a signal processing flow graph.

In your case, it'd look like this:

Vector Source -+> high pass filter --> threshold --> [custom block]

which  you can click together graphically in GNU Radio companion and let it generate Python, or you can directly write in Python.

The [custom block] is from the top of my head something you'd have to write yourself. You can do that in Python or C++, but Python will be slower to run. I'd still do it in Python for a first prototype:

When the first input is 0, it simply consumes the the input (effectively dropping that data). When it's 1, save the number of the first 1-sample in a row. Then wait until the 1s stop and safe locally something

self._list_of_signals.append( (first_item_number, last_item_number ) )

I think you get the idea.

You feed in the data into the vector source in your python framework, e.g. with the data from your database, and then let the flow graph run, i.e. process all the data it has.

When it's done running, you ask your custom block to give you the list of timestamps. With their first sample number times the sampling rate, you get a time stamp, and then you just put (start, stop) into your database of "trimmed" data.


How to identify stable deltas in a noisy signal

This sounds like the "change detection" problem of Basseville and Nikiforov, with the parameter after the change unknown (see section 2.4) of the linked PDF.

See, for instance, example 2.4.3:

The usual algorithm for dealing with this, and the one analyzed in the above-linked book is the CUSUM (Cumulative Sum) algorithm.


Plane Settings of the Matched $z$-transform Method

The $s$-plane is the complex plane associated with the Laplace transform, i.e., with transfer functions of continuous-time systems, whereas the $z$-plane is the complex plane associated with the $\mathcal{Z}$-transform, i.e., with transfer functions of discrete-time systems. The matched $Z$-transform is one way of transforming continuous-time systems to discrete-time systems.

Take a look at this answer for more information about transformations from continuous-time to discrete-time.


The proof of associativity of discrete convolution relies on the assumption that multiple infinite sums can be evaluated in any order. This is not true if some of the involved sequences do not converge absolutely, which is the case for the given sequences $x_1[n]$ and $x_2[n]$. Note that the convolution sum $x_1\star x_2$ does not converge, i.e., $x_3\star (x_1\star x_2)$ gives yet another (infinite) result.

In continuous time you have the same problem. Associativity of continuous convolution relies on Fubini's theorem for switching the order of integration. If the assumptions of Fubini's theorem are not satisfied for the given functions, convolution is not associative.


Tricks for plotting the magnitude of a DTFT?

So the example that you've listed appears to be somewhat odd. Assuming Ts is the time step, you're saying that the continuous signal x(t) is dependent on the time step you've chosen. That just seems odd.

It also seems that you've formatted your question in the sense that you understand the Fourier Transform (FT), but you're just having a difficult time relating the FT to the associated Discrete Time Fourier Transform (DTFT) after sampling the continuous signal. I'm going to speak at a relatively high level to explain how to go from the FT to the DTFT.

The first thing you should realize is that the DTFT will be related to the FT of the continuous signal. But how so? The sampling theorem (which is the most fundamental concept in DSP) tells us that the bandwidth of a discrete signal (i.e. sampled signal) is equal to that of the sample rate$^1$. We typically write this as -Fs/2 to Fs/2 where Fs is your sample rate. Furthermore, we know that analog frequencies that are integer multiples of the sample rate will be interpreted as equivalent frequencies. This is aliasing.

For example, say Fs = 1kHz. Let's randomly select a frequency of 0.5kHz. All frequencies that are integer multiples of the sample rate at this frequency will be seen as 'the same thing' when sampled. 0.5kHz, 1.5kHz, 2.5kHz, etc. etc. will all appear to be 'the same signal' when sampled. This is also true of other frequencies such as 0.25kHz. In this case, 0.25kHz, 1.25kHz, 2.25kHz, etc. will all appear to be the same. Note that this phenomenon happens at EVERY frequency.

This is what we refer to as aliasing: different frequencies be interpreted as the same signal due to sampling. A common misperception is that aliasing only happens in one direction: that being that high frequencies are aliased onto lower frequencies. However this isn't true. Looking at our above examples, you can see that aliasing only comments on the fact that different signals are interpreted the same. It doesn't say anything about direction. Thus you could also state that lower frequency signals are aliased onto high frequency signals. Because of this aliasing of high frequency signals onto low frequency signals, and low frequency signals onto high frequency signals, the frequency spectrum of sampled signal becomes periodic. So let's put this together into a simple list of steps on how to get the DTFT from the FT.

1) Obtain a copy of the FT of the continuous signal. You will need to do this analytically.

2) Slice up the FT along the frequency spectrum at intervals of Fs. Let's say your first interval would range from -Fs/2 to Fs/2 (this is centered at zero).

3) After you've 'cut up' your FT into intervals, you can take all of those and add them together at a single interval (i.e. do the super position of them). By doing this step, we are essentially introducing aliasing into one interval.

4) Now repeat step 3 for every interval. By doing this, you're accounting for the fact that aliasing doesn't just occur from high frequency to low frequency, but it essentially occurs everywhere!

5) Look at what you've done. By doing the super position at every interval, you'll see that your frequency spectrum is now periodic (and still continuous!).

This is essentially what the DTFT is doing. Visualizing these steps is very useful, I might add illustrations if I get time. Now assuming that you've got a plot of the FT, you now have a plot of the DTFT.

[1] It would be more accurate to say that the relevant bandwidth is that of the sample rate.


Distortion of OFDM signal when using RRC filter

I checked your code in my PC, you need just to delete the delay added before the filter. For example, you can use:
U_aft_fil = U_aft_fil(fil_delay+1:end);

Then when filtering it again at the receiving side, you delete it again :
U_r_fil = U_r_fil(fil_delay+1:end);

Good luck


Does downsampling increase the resolution of frequencies?

Downsampling by a factor of $N$ in time-domain means that you throw away $N-1$ samples from $x[n]$ for every $N$ samples. In frequency domain this creates $N$ shifted copies of the original spectrum and expansion of frequency axis. The shifted copies are shifted by $\omega = 2\pi \frac{k}{N}, \ k = 0,1,2,...,N-1$. So, the DTFT of downsampled sequence $x_D[n]$ is basically given by:
$$X_D(e^{j\omega}) = \frac{1}{N}\sum^{N-1}_{k=0}X\left(e^{j\left(\frac{\omega}{N} - 2\pi \frac{k}{N}\right)}\right)$$
This is why downsampling can cause aliasing and to avoid aliasing we need to Low pass filter the original sequence to suppress sprectrum outside of $|\omega| > \frac{\pi}{N}$.
Now, if the original signal was appropriately oversampled, then the DTFT would have been occupying only $\omega \in [-\frac{\pi}{N}, \frac{\pi}{N}]$. And, the downsampling will expand the spectrum by a factor of $N$ without aliasing and it will give an impression that the frequency resolution has increased. Because you are now seeing the same spectrum which extended from $\omega \in [-\frac{\pi}{N}, \frac{\pi}{N}]$ stretched to $\omega \in [-\pi, \pi]$.


It's probably more straightforward to first compute the total transfer function in terms of $H(z)$, then plug in the actual expression for $H(z)$, and from that write down the difference equation:
$$Y(z)=X(z)+\alpha z^{-M}H(z)Y(z)\tag{1}$$
From $(1)$ we can derive the total transfer function:
$$\frac{Y(z)}{X(z)}=\frac{1}{1-\alpha z^{-M}H(z)}\tag{2}$$
Plugging $H(z)=(1-\lambda)/(1-\lambda z^{-1})$ into $(2)$ gives
$$\begin{align}\frac{Y(z)}{X(z)}&=\frac{1}{1-\alpha z^{-M}\frac{1-\lambda}{1-\lambda z^{-1}}}\\&=\frac{1-\lambda z^{-1}}{1-\lambda z^{-1}-\alpha(1-\lambda)z^{-M}}\tag{3}\end{align}$$
Eq. $(3)$ is equivalent to
$$Y(z)=\big(1-\lambda z^{-1}\big)X(z)+\big(\lambda z^{-1}+\alpha(1-\lambda)z^{-M}\big)Y(z)\tag{4}$$
from which you can directly see the corresponding difference equation.


Looking for the Concept About All In One Curve Fitting

You can always augment the matrices to do so.

Let's assume the first model is given by:

$$ {y}_{1} = {H}_{1} * {\theta}_{1} $$

The second model is given by:

$$ {y}_{2} = {H}_{2} * {\theta}_{2} $$

The third model is given by:

$$ {y}_{3} = {H}_{3} * {\theta}_{3} $$

If we assume the number of parameters of the model are the same, namely $ {\theta}_{1}, {\theta}_{2}, {\theta}_{3} \in {\mathbb{R}}^{d} $ and $ {y}_{i} \in {\mathbb{R}}^{{n}_{i}} $ by the following augmentation:

$$ y = H \theta $$

Where $ y = {\left[ {{y}_{1}}^{T}, {{y}_{2}}^{T}, {{y}_{3}}^{T} \right]}^{T} $, $ H = {\left[ {{H}_{1}}^{T} \mid {{H}_{2}}^{T} \mid {{H}_{3}}^{T} \right]}^{T} $ and $ y = {\left[ {{\theta}_{1}}^{T}, {{\theta}_{2}}^{T}, {{\theta}_{3}}^{T} \right]}^{T} $ you can get:

$$ \hat{\theta} = {\left( {H}^{T} H \right)}^{-1} {H}^{T} y $$

Extracting the the corresponding parameters is easy.

If the parameters vectors are not the same dimension, you can use zeros to make everything work.


In a real system you actually need two DACs: one for the real part and one for the imaginary part. The implementation of the DAC model depends on what should be included in your simulation. There are two main aspects: quantization and pulse shaping.

Quantization - if quantization effects of the DAC should be modelled, the floating point output of IFFT/CP insertion has to be mapped to a discrete number of values. Quantization can, for example, be implemented by scaling and rounding. Say your DAC has a vertical resolution of 8 bit. Then the IFFT output signal should be scaled such, that its maximum value is 127. Afterwards the signal is rounded to integer numbers. Note that OFDM signals are often clipped before quantization to reduce quantization noise, especially for lower resolutions (1).
Pulse shaping - if a physical channel should be modeled (for example a wireless channel), the analog output of the DAC must be simulated. This is implemented by converting every input sample into an impulse that is weighted with the sample's value. In the simulation the impulse itself consists of $L$ samples, of course, that mimic the analog signal. $L$ is the oversampling factor (not to be confused with the OFDM oversampling that is achieved by zero padded subcarriers). The simplest impulse shape is a rectangle, i.e. every sample in the simulation is just repeated $L$ times.

If the channel model is very simple, for example an AWGN channel, DAC and ADC can be omitted. It really depends on how realistic the simulation model should be.

(1) Michael Bernhard et al.: Analytical and numerical studies of quantization effects in coherent optical OFDM transmission with 100 Gbit/s and beyond


Is it possible to use a control scheme with a delay as prefilter?

1 - It can be a good idea to put a filter $C_r(s)$ in order to smooth the reference. For example, if you have a speed regulator, and the reference jumps from 100 km/h to 200 km/h, having a low-pass filter will smooth the actual reference fed to the closed-loop. It could prevent excitation of unmodeled dynamics, actuator wind-up, etc. That being said if $C_r(s)$ is a simple delay, it has no purposes as far as I can see.

2 - $C_r(s)$ only affects the reference fed to the closed-loop. It will not affect stability of the closed-loop like you mentioned. It also has no effect on disturbance rejection or errors due to modelization (such as unmodeled dynamics).


Impulse and step response of a system that contains cascaded second order sections

Not sure I follow completely, but is this the question: how to obtain the impulse and step response of a filter that is realized as a cascade of second or sections? or are you asking for working code that implements this?

Your 8th order filter is decomposed into 4 second order sections. The ordering is not really important here. For any input $x$ to your 8th order filter you compute the output $y$ by first filtering the input with the first second order section $SOS_1$ to get output $s_1$, that is $s_1 = SOS_1*x$ (the $*$ means processing the $x$ with a second order difference equation with coefs from $SOS$). Then that output is made input to the second SOS, that is $s_2 = SOS_2*s_1$. So now you have $s_2 = SOS_2*(SOS_1*x)$. And then continue this pattern (output from one section is input to the next) through all sections.

If you understand this, then it should also be rather obvious that you obtain the impulse response or the step response by letting an impulse or a step waveform be input to your filter and the output $y$ contains the response.

I have not looked at your code. But try to start with a simpler case first just to verify your code and your understanding of the concept.


Why we take different mfcc components in different programs?

The number of coefficients is just one of the many parameters to tweak when designing a classifier. For speech analysis it turned out that usually 8-14 coefficients are sufficient, and increasing the number of coefficients beyond 14 does not improve the performance for most applications. The exact number you choose will depend on your task at hand. Most systems I've worked with used 12 coefficients, and that was for speech recognition. If you're not sure, you can simply run some experiments and change the number of MFCCs in a reasonable range to see the effect it has for your application. But I have no doubt that you'll end up in the range that I've mentioned above.


Why the error probability is given by this expression?

A whole chapter of a book can be written to answer the "Why" in this question but the
short answer is that if the signals used are an antipodal signal set (meaning
bits $0$ and $1$ are transmitted using signals $s_0(t)$ and $s_1(t)$ respectively
where $s_1(t) = -s_0(t)$ and the signals are of duration $T$ and energy $E$,
then with matched filtering, the receiver makes a decision that $0$ was transmitted
or a $1$ was transmitted according as the sample value from the matched filter
is positive or negative.  The sample value is a random variable $X$ whose
distribution is $\mathcal N\left(\pm \sqrt{E}, N_0/2\right)$ depending on whether
a $0$ or a $1$ was transmitted. If $X$ has mean $-\sqrt{E}$ when a $1$ is transmitted,
then the probability of error is
$$\begin{align*}
P_{e,1} &= P\{X > 0 \mid 1~\text{transmitted}\}\\
&= P\left\{\mathcal N\left(-\sqrt{E}, N_0/2\right) > 0\right\}\\
&= Q\left(\frac{0 - \left(-\sqrt{E}\right)}{\sqrt{N_0/2}}\right)\\
&= Q\left(\sqrt{\frac{2E}{N_0}}\right).
\end{align*}$$
Similarly, $X$ has mean $\sqrt{E}$ when a $0$ is transmitted,
and so the probability of error is
$$\begin{align*}
P_{e,0} &= P\{X < 0 \mid 0~\text{transmitted}\}\\
&= P\left\{\mathcal N\left(\sqrt{E}, N_0/2\right) < 0\right\}\\
&= \Phi\left(\frac{0 -\sqrt{E}}{\sqrt{N_0/2}}\right)\\
&= \Phi\left(-\sqrt{\frac{2E}{N_0}}\right)\\
&= Q\left(\sqrt{\frac{2E}{N_0}}\right).
\end{align*}$$
Finally, using the law of total probability, the probability of error is
$$P_e = P_{e,0}P\{0~\text{transmitted}\}
+ P_{e,1}P\{1~\text{transmitted}\} = Q\left(\sqrt{\frac{2E}{N_0}}\right)$$
regardless of the probabilities of transmitting $0$'s and $1$'s.
For more details such as why the matched filters outputs
are as stated above, see, for example, this lecture
note and
this one
of mine.


Determining if a system involving integrals is time varying or time-invariant

The fixed bound at 0 in $\int_0^t$ indeed suggests a time-variance.You only need a counter-example for $S_1$, like the constant $x(t) = 1$.  $y(t+1) = t+1$, while $S_1(x(t+1)) = S_1(x(t)) = y(t)$.

For $S_2$, the two bounds moving at the same time are promising, and a standard variable change shows the time invariance.


When discussing the z-transform, it is important to note that there are two different frequencies that are innately involved. In many text books, the two frequencies are denoted by $\omega$ (discrete-time) and $\Omega$ (continuous-time). We'll use this convention.
The z-transform uses the "discrete frequency" $\omega$ which is derived from the "analog frequency" $\Omega$ via sampling at a rate $T_s$ as
$$\omega = {\Omega}T_s$$
If we're on the unit circle, which your question implies as the starting condition, then the $\sigma$ term in $s = \sigma + j\Omega$ goes to zero and so
$$s = \sigma + j\Omega \Rightarrow s = j\Omega$$
We can now relate the z-transform as the discrete-time equivalent of the Laplace transform and so
$$z = e^{sT_s} = e^{j{\Omega}T_s} = e^{j\omega}$$
Which is exactly the substitution made between the Fourier transform of a discrete sequence $x[n]$ and it's z-transform:
$$X(e^{j\omega}) = \sum_{-\infty}^{\infty}x[n]e^{-j\omega n}$$
$$X(z) = \sum_{-\infty}^{\infty}x[n]z^{-n}$$
I have never seen $z$ be defined by $z = e^{s}$, which implies the use of the continuous-time frequency and is not how the z-transform is defined.
So the confusion might come down to abuse of notation between how $\omega$ and $\Omega$ are defined in the sources you are referencing.


Practical wideband digital beamforming for large arrays in radar applications

I've not worked on the design of such systems before, but I think your notions are on the money. Specifically, yes, beamforming arrays do have RF front ends that are replicated many times. The complexity of contemporary phased array radars is astounding in this regard; there are designs that have hundreds of individual antenna elements in them with impressive levels of control of the array response using various signal processing techniques.

And as you suspected, yes, this sort of approach is not cheap. Gigasample-class ADCs are available commercially in the few-thousands-of-dollars range, but it's possible that the custom, low-quantity RF front ends used in systems like this would dwarf that cost. Even so, radars with this sort of capability are often found as subsystems in very expensive larger systems (like a multi-hundred-million-dollar fighter jet).

As far as the backend digital signal processing goes, that's a pretty mature market that has developed over the past few decades. The main goal is processing density: getting the maximum number of FLOPS into the smallest volume. After all, such radars are often used in space-constrained applications like aircraft. Therefore, you'll see a lot of the processing done on custom FPGAs and/or single-board computers that can be stacked compactly into standardized chassis assemblies (like VPX or CompactPCI).


Continuous time-domain electrical signal that only has a finite number of possible values

Well, a DC signal for sure.

Other than that, there is probably more context from the book required (especially about the assumptions).


Testing if a random signal was modified or if it contains noise

No.  You have to start with some prior knowledge of the signal and of the noise.  Your signal may be my noise, and visa-versa.


Generate loopable sine wave cycle for given frequencies

Your clicks are coming from two sources. The wavetable and the hardware.

If you are looking to create a simple little instrument that produces tones by playing back a wavetable, then there is no way to fix the size of the "playback window" to reproduce tones without clicks.

The frequencies that compose the chromatic scale follow a geometric progression that achieves doubling of the frequency in 12 steps. The frequency for each tone is produced by multiplying the previous tone by something like 1.05. That factor is worked out by: Let's take $A_F = 440Hz$ and its octave at $A_O = 880Hz$. Each one of the tones is $A_{n+1} = A_n * r$, so $A_O = A_F \cdot r^{12}$. But we know that $A_O = A_F \cdot 2$. Therefore, $r^{12} = 2$ and $r \approx 1.059463094359$

So, you might fit $A_F,A_O$ neatly in a wavetable, but $A_F \sharp$ now needs $\approx 1.0594$ more samples to fit in neatly. Not only there is no integer there but even if you could reproduce the tones in neat windows you would get problems with the tempo, the duration of each note in a musical piece, because $A_F$ is shorter than $A_F \sharp$.

The "classic" way around this is to setup a high resolution wavetable of a single sinusoid which is then read back at "different speeds", with interpolation. So, $A_F$ is read back at $r^0 = 1$ speed (and practically no interpolation) but $A_F \sharp$ is read back at $r^{1} \approx 1.0594$ speed and a bit of interpolation when the "needle" has to land between known samples.

BUT!, this has nothing to do with your output stream. In other words, you have to set-up (or adopt) a stream system where samples are continuously being pushed to the sound card and you have a "global" sense of time.

Once upon a time, you did this by setting up three pointers (that's in C, C++) to three buffers. Let's call them $H,D,L$. You set up the sound card to record in $H$ and playback $L$ and while that was going on, your main code was processing buffer $D$. When the sound card called its callback to signify that the recording ended, you would "circularly switch" the buffers around so that they now pointed to $L \rightarrow D, D \rightarrow H, H \rightarrow silence$ and re-set up the sound card to do the same job while you were now processing the $D$ which is now your freshly recorded $H$. Those buffers had a length, say for instance 1000 samples and the task then was becoming to find a sweet spot that balanced length of buffer with responsiveness.

But as you can see, in this system, samples are being pushed out to the sound card continuously in blocks of $bufferLength$ samples and you don't really care if what is being played back fits the buffer. All that you do is iterate your models (e.g. the wavetable), produced a buffer of samples, freeze the models, playback that buffer and move again from the top.

Now, obviously, the above includes real time input from audio sources too. If all you are doing is a tiny little synth, you could have two buffers and switch around between them, working on one while the other is being played back.

But even if you did this with some off the shelf cheap sound card in the market, you still got clicks although, you have worked so hard to perfectly align the zeroes between blocks. This is because of something called latency.

In other words, just because your program called some Operating System Application Programming Interface function that said "Record NOW!" (or "Play NOW!), it doesn't mean that this command will propagate through the driver and eventually the audio hardware and recording will start NOW. Recording (and playback) will commence some $NOW! + latency$ units of time later and in off the shelf sound cards, this used to be an awful time later. So awful that you could hear the silence between calls.

So, the solution there was to move over to something like ASIO which eventually packaged this whole workflow into a complete solution all the way from the software calls to the hardware.

As an example, an 8-track external sound card can do 48kHz Stereo playback / rec over buffers of 256 samples (with absolutely no clicks). This means, that sound is processed in blocks that last approximately 5ms. The human ear begins resolving sounds as distinct when they are apart for more than 40ms, so this 5ms is really safe. It means that we can plug an electric instrument in the sound card and process it, in real time like it went through an effects box and the performer will not be noticing any difference between what they play and what they hear.

If you tried to make $H,D,L$ the size of 256 and try to do the same with an off the shelf (relatively popular) sound card, first of all you would g/et/ /an/awf/ul/am//ount/of/cli/k/s// (one every 5ms) and secondly, the driver and operating system would be so much overwhelmed by the amount of function calls that soon you would be greeted with Seg Faults, blue screens, system freeze, reprepreprepreprepreprepreprepreprereprepreprepreprepreprepeated waves out of synsynsynsynsynsynsynsynsynsynreprepresync and other amusing side effects.

So, bottom line, keep working on the DSP side of your synth but make sure that your sound hardware is up to the task as well if you really want to get rid of the clicks.

Hope this helps.


That sounds a lot like you shouldn't be using the DFT, but rather a parametric spectrum estimator.

For example, there's spectrum estimators that actually give you a function to evaluate at real-valued points rather than a vector of a sampled spectrum (your sentence "DFT, unquantized" doesn't make all that sense to me – the DFT is discrete, inherently, so what you describe as search algorithm only works if you do a massive DFT and then start somewhere in the and look for peaks from these points on).

My go-to example here is either

ESPRIT, which, given knowledge about the number of significant signals will just give you a vector of real-valued frequencies of peaks or
ROOT-MUSIC (dunno if that's the official term, it's what we used, based on [1]) which uses the pseudospectrum function generated by the MUSIC algorithm and finds the signals by looking for roots of a specifically crafted polynomial [1].

[1]: Barabell, Arthur J.: Improving the resolution performance of eigenstructure-based direction-finding algorithms. In: Acoustics, Speech, and Signal Processing, IEEE International Conference on ICASSP ’83., Band 8, April 1983, available online.


Decoding Matrix Encoded Surround Sound Through Convolution

Suppose that we had two channels encoded in some form in two signals according to a matrix.  Recovering the original signals would be straightforward using the matrix inverse (or impossible, if the inverse did not exist).  The situation is no different here, except that the matrix is not square, so you cannot use the conventional matrix inverse.  You can, however, use a pseudoinverse or generalized inverse to reconstruct the original signals.

Other matrices which are not equal to an inverse may work also succeed, and may in fact sound better.  I don't know the psychoacoustic situation, but it may more sense to do some transform based on the position of speakers for example. Certainly, the inverse seems like a good place to start.

I rewrote your encoder table as a 2x5 input matrix $ E $, and used Wolfram Alpha to determine the inverse, $ D = E^{-1} $.  Here is the 5x2 result it provided, using $ i $ as the imaginary unit, with $ j = i $ and $ k = -i $.

$$
D = \frac{1}{18(46+6\sqrt{2})} \begin{bmatrix}
405 & -3(27 - 36\sqrt{2}) \\
-3(27 - 36\sqrt{2}) & 405 \\
6(18+27\sqrt{2}) & 6(18+27\sqrt{2}) \\
-i(27\sqrt{3} + 99\sqrt{6}) & i(63\sqrt{3} + 27\sqrt{6}) \\
-i(63\sqrt{3} + 27\sqrt{6}) & i(27\sqrt{3} + 99\sqrt{6})
\end{bmatrix}
$$

As one would expect, the decoding matrix shows left/right symmetry.

Now we can express the encoding/decoding operations as:
$$
\begin{bmatrix}
L_{encoded} \\
R_{encoded}
\end{bmatrix} = E \space{} \begin{bmatrix}
L \\
R \\
C \\
L_{rear} \\
R_{rear} \\
\end{bmatrix}
$$
$$
\begin{bmatrix}
L' \\
R' \\
C' \\
L_{rear}' \\
R_{rear}' \\
\end{bmatrix} = D \space{} \begin{bmatrix}
L_{encoded} \\
R_{encoded}
\end{bmatrix}
$$

Now you can easily obtain formulas to express the three front components as functions of the encoded left and right signals.  For the two rear components, you have an equation in the time domain with a complex component.

Looking at the left rear component, and converting constants from the matrix $ D $ to numeric form:
$$
L_{rear}'[n] = -i0.2949 L_{encoded}[n] + i0.1787 R_{encoded}[n]
$$

The Hilbert transform can be used to remove the imaginary unit, because this has the same effect on the frequency domain of performing a ±90º phase shift.  The Hilbert transform only needs to be applied once, since convolution ($ * $) is distributive and commutative and all other good things.
$$
L_{rear}'[n] = (-0.2949 L_{encoded}[n] + 0.1787 R_{encoded}[n]) * H[n]
$$

The operating principle, then, is that the encoding technique functions in part by encoding the left/right rear components in quadrature.

The Hilbert transform is infinitely long, so it must be truncated and expressed as either a FIR filter, or depending on your circumstances, a FFT-IFFT operation.  A causal Hilbert transform will delay the rear channels, so your other filters will need to be delayed by some amount to compensate.

For example, if $ k $ is the group delay of a causal approximation of a Hilbert filter $ H_k[n] $ then you delay must delay the three front channels by $ k $ as well:
$$
C'[n] = (0.3437L_{encoded}[n] + 0.3437R_{encoded}[n]) * \delta[n-k]
$$

In practice, all of the output signals depend on the two input signals.  That is, you take two input channels and produce one output channel five times.  Your JACK pipeline will need to accommodate this.


FIR filter design with flat passband but equi-ripple stop-band

I suggest you the Parks-McClellan method. In Matlab you can use $\tt{firpm}$.

Matlab's (now obsolete function) $\tt{remez}$ also uses this scheme.

The FIR filter is optimally designed to approximate e.g. a Chebyshev filter such that the maximum error between the filter's response and the desired response is minimized.

For more information see this and the examples therein.


Separate speech signals that are bandlimitted

As mentioned in comments that the 4Khz and 8 Khz bandlimited files are kept separately but in same folder and problem is to segregate them into separate folders.

One can use scripts written in Matlab or Python for doing this, procedure could be as follows :

Read each file (this could be easily automated in matlab/python) from folder.
split each file into frames, take FFT of frames.
for the 4 Khz band limited files the magnitude of upper half of FFT coefficients will be very less, but for 8Khz magnitude of upper half of FFT coefficients will be significantly high.
the above logic could be automated by keeping fixed threshold on magnitudes of FFT coefficients.
based on above threshold, files could be moved to separate folders.


Why is an interleaver often placed between an outer block code and an inner convolution code?

As the other answers mention, the use of FEC results in post-decoding errors occurring
in bursts. Indeed, this happens regardless of whether the code is a convolutional code
or a block code.  With a $(n,k)$ block code, the decoder output ($k$ bits)
from the decoding of one
received word is (hopefully with high probability) completely correct, or it has
an unknown number of errors in it that can be regarded as a burst error of length $k$.
With a convolutional code, the decoder output is mostly correct as the decoder
finds the correct path through the trellis, but occasionally the decoder's
chosen path deviates from the correct path and later rejoins the correct
path during which time there is a burst of errors. In contrast to
block codes, there can be multiple
isolated burst errors (of variable lengths) in a single transmission
using a convolutional code. Also, the number of data bits in a single
transmission is far larger than the typical values of $k$ for a
block code.

The idea behind a concatenated
coding scheme is that with an inner code suited to the physical channel,
and an outer code over a very large symbol alphabet, we can make
the burst errors in the inner decoder output look like a single
symbol error to the outer code.  This is important because the
outer code should be a very high rate code because the net rate
is product of the inner code rate (more or less determined by
the channel and the link budget) and the outer code rate.
Unfortunately, outer codes over very large alphabets are very
expensive to implement, and so interleaved Reed-Solomon codes
over smaller alphabets (often GF$(2^8)$ ) are used (with
interleaving at the symbol
level, as Jim Clay points out). Because of the interleaving, the burst
errors in the inner decoder output become single symbol (byte) errors in the
received words of the interleaved Reed-Solomon code.

All the above is mostly a rehash of what the answers by Bryan and
Jim Clay have already said, but I wish to point out the following.

  Interleaved Reed-Solomon codewords can be decoded much more
  efficiently and with smaller delay if they are not de-interleaved
  first.

A Reed-Solomon decoder that can decode interleaved codewords is different from
the off-the-shelf standard Reed-Solomon
decoders that are available, and the use of such a decoder might not be
feasible if the development team does not have control of this aspect of
the design. But, if such a decoder is used, the de-interleaver
can be moved from the between the inner decoder and outer decoder to
just after the outer decoder. The de-interleaver is also smaller since
it has to deinterleave a $K\times L$ array instead of a $N \times L$
array (for a $(N.K)$ Reed-Solomon code interleaved to depth $L$.

  If a delay-scaled Reed-Solomon encoder is used along with the
  Reed-Solomon decoder for interleaved codewords described above, the interleaver
  at the transmitter and the de-interleaver at the receiver
  can be eliminated entirely.

The output of the delay-scaled
encoder (see also
this paper
which is unfortunately behind a paywall)
is a set of interleaved Reed-Solomon codewords, but
is not the same sequence of bytes that one would get from
doing a standard Reed-Solomon encoding followed by interleaving.
So, no further interleaving is necessary.
The output of the interleaved Reed-Solomon decoder is the
same byte sequence in the same order that went into the delay-scaled
Reed-Solomon encoder, and so no de-interleaving is necessary
at the decoder, either.


Code generation for `sum` function in MATLAB

Because the size of y1(j1:k1) and y2(j2:k2) must be the same then:
$$
k_1-j_1 = k_2-j_2
$$
and you will also have something like: $j_1+N=j_2$ then you could write a for loop like this:

xc(i) = 0
for n=j1:k1
  xc(i) = xc(i) + y1(n)*y2(n+N)
end


PC-to-PC data communication over Sound FSK / PSK / DSSS

This is actually a really tough problem because of the channel characteristics. Most computer speakers have fairly limited bandwidth, have significant non-linearities and the room acoustics are often time variant.

Life becomes A LOT easier if you can just run a cable from the headphone output of one PC into the line input of the other.


I think I've found the answer. Please correct me if I'm wrong. First of all, I've made a silly mistake  $$\int_{-\infty}^{+\infty}|y(t)|dt = \int_{-\infty}^{+\infty}|Y(s)|ds$$which is clearly false. Let $y(t) = x(t)\star h(t)$. We have $$E_y = \int_{-\infty}^{+\infty}|y(t)|^2dt = \int_{-\infty}^{+\infty}|Y(s)|^2ds = \int_{-\infty}^{+\infty}|H(s)X(s)|^2ds = \int_{-\infty}^{+\infty}|H(s)|^2|X(s)|^2ds$$Also we have $$|H(s)| = \left|\int_{-\infty}^{+\infty}e^{-2\pi ist}h(t)dt \right | \le \int_{-\infty}^{+\infty}|e^{-2\pi ist}h(t)|dt =  \int_{-\infty}^{+\infty}|h(t)|dt \lt\infty$$So  $\exists M\in\mathbb{R}:\ \ |H(s)|\le M$ for all $s$. This means that $|H(s)|^2\le M^2$ and then $$\int_{-\infty}^{+\infty}|H(s)|^2|X(s)|^2ds\le M^2\int_{-\infty}^{+\infty}|X(s)|^2ds$$By assumption $$E_x = \int_{-\infty}^{+\infty}|x(t)|^2dt=\int_{-\infty}^{+\infty}|X(s)|^2ds$$The result is $$E_y \le M^2 E_g$$


Meaning of Rect and Train of Rect Spectra

Why does repeating periodically a signal means removing some frequency intervals from the single pulse spectrum?

Because your signal is periodic. A periodic signal only contains frequencies that are integer multiples of the signal's repetition frequency. In other words, if the signal is periodic all of it's constituent parts must be periodic as well.

In time domain, frequency may be seen as the speed of variation of the signal.

Not really. "speed of variation" would be the first derivate of the signal. It's somewhat related to higher frequencies since the Fourier Transform of a differentiator is proportional to the frequency itself but it's not the same thing. If you look at a square wave you also need high frequency to keep the flat parts nice and flat.

So, my question is a about the physical meaning of this.

That's always a tricky proposition. One possible answer is "there is no physical meaning since it's physically impossible to have a signal that's infinite in time". Sine waves are a useful mathematical concept but they DO NOT EXIST in reality. Of course, for practical purposes you can always generate something that's "close enough" for your specific application.
Another explanation
Let's look at what happens if you repeat the rect just once. In the time domain you can represent this as a convolution with two pulses, i.e.
$$h(t) = \delta(t) + \delta(t-T_0) $$
So you add a delayed copy to the original signal. In the frequency domain this looks like
$$H(\omega) = 1 + e^{-j\omega T_0}$$
The delay turns into a phase shift so you add add a phase shifted copy of the spectrum the original one. If the phase shift at a certain frequency is 0 degree, they will just add. But if the phase shift at a different frequency is 180 degrees they will actually cancel. So some frequencies will be amplified others will be cancelled that's called a "comb filter" a graph of the amplitude with these notches looks like the teeth of a comb.
The more you repeat the original rectangle, the more "combing" you get and eventually all frequencies are gone except the ones that line up exactly with the period of the repetition.


DFT Frequency domain analysis and interpolation

A while back, I remember some professor saying that for the N point DFT Frequency Domain, the values for k around N/2 (halfway point for the spectral domain - those around π) are considered high frequencies, and the end points (0, N-1) (those around 0, 2π) are considered low frequencies. I am having trouble seeing this mathematically given the definition for computing the coefficients of the DFT.

You can define the DFT however you like, as long as in effect the same frequencies are covered - remember, $f(x) = e^{jx}$ is $2\pi$-periodic, so that, for example, the sum over $k=0,\ldots, N-1$ is the same as the sum over $k=-\left\lfloor\frac N2\right\rfloor,\ldots,0,\ldots,+\left\lfloor\frac {N-1}2\right\rfloor$.

  The examples I have seen seem to decide to pad the coefficients around the middle, i.e. for N=4 and a desired sampling of 2N, [1 2 3 4] become [1 2 0 0 0 0 3 4]. Is there a reason why the middle is targeted and not the end like [1 2 3 4 0 0 0 0]?

Again, different definitions of the DFT lead to different forms for zero-padding here.


Performance metrics for analysing the performance of various filters

They are getting values like that because they may be reading the image as double while calculating the metrics, you can use im2uint8 to convert the double into uint8 and estimate the metrics. you will get the desired output. Though it does not mean that their values are incorrect.


Your math is correct, for a delay of 2 seconds, the first null in the PSD should be at 0.25 Hz.

I'm guessing the problem is in the code you use to generate the signals and/or the PSD. The sample rates seem awfully low. At a sample rate of 0.25 Hz, your 2 second delay is just half a sample, that means you would have to implement a fractional delay to generate the signal. T


Quadrature subsampling result not as expected

Your calculation is correct; sif has a harmonic at 308 kHz. Given the way you have set up the problem, you can't get that exact value using the FFT because you don't have a frequency bin at that frequency. According to my calculations, the nearest bin is at 305 kHz, and I do see a peak at that frequency. Since there is no bin at the signal's frequency, the energy will leak: bins in the neighborhood of the peak will also get some of the energy of the 308 kHz signal.

Leakage is very hard to avoid. Assuming a fixed sampling frequency, you can calculate the number of samples you need to have a bin at exactly 308 kHz. However, it's very likely that you will then have a non-integer number of cycles of sif, so you will have leakage again. This effect can be reduced by using windowing.

There are many related questions in this website (for example, this one). Search for 'fft bin leakage' to find many more.


Remove Checkered Noise (Pattern) from an Image

You could model the FFT probabilistically and use your prior knowledge about the FFTs of natural images, such as the one of the left, to model a prior distribution with which to form a MAP estimate.


Fourier transform as the integral of a parameter multiplied by an homogeneous wave

given your definition, and the definition of the Fourier transfer:
yes. That is literally the definition:
$$\int_{\mathbb R} s(x) e^{i2\pi f x} \mathrm dx$$
By comparison, with $k= 2\pi f$, you get the Fourier transform. Every book uses that definition, so you'll find this "interpretation" in every book that introduces the Fourier transform. (This really feels like exactly the same $k$-space / location transform I learned in my first solid state physics lecture, so I assume every modern solid state course for physicists does pretty much that, but more in-depth.)


Each point in surface has a normal represented by 3x1 vector so transpose of it is 1x3 . Considering radiance from the source is represented by normal to the source it will be represented by a 3x1 vector.


I've used the Liu-Liang-Weisberg paper before, and it worked pretty well. There's some example code on this page.

In short, they propose that squaring wavelet coefficients and dividing them by their respective scales normalizes the spectrum. The paper provides proof and motivation.


First, your period is $T = 1/25s=0.04s$ which means you do not need to take integral of 3 periods for getting fourier series. A plot of the DFT is given below when I took only first $T=0.04s$ is given by

As you can see, all even multiples of $25\text{Hz}$ is zero, while odd multiples of $25\text{Hz}$ is decreasing with increasing frequency. From the plot $|X[5]|^2 = 0.002914$, while $|X[7]|^2=0.0008386$. Even multiples of $25\text{Hz}$ have very low values ($\lt 10^{-12}$). Taking DFT is not the correct method, so read on.

The above values are expected because fourier series of your even symmetric signal with non-zero DC has all $b_k =0$. For computing $a_k$, your signal $x(t)=2t/T$ for $0 \le t \le T/2$ and $x(t)=-2t/T-2$ for $T/2 \le t \le T$. I am not going into details of how to compute fourier series coefficients but I will give hints
$$
a_0=\frac{2}{T}\int_0^T x(t)dt\\
a_k=\frac{2}{T}\int_0^T x(t)\cos(2\pi kt/T)dt \\
=\frac{2}{T}\int_0^{T/2} 2t/T\cos(2\pi kt/T)dt + \frac{2}{T}\int_{T/2}^T (-2t/T-2)\cos(2\pi kt/T)dt\\
=\frac{2((-1)^k-1)}{\pi^2k^2}\\
a_k=\frac{-4}{\pi^2k^2}\text{  for k odd}\\
a_k=0\text{ for k even}
$$
So as you can see, the even harmonics are all zeros, while odd harmonics decrease by a factor of $1/k^2$. So
$$a_5=\frac{-4}{\pi^2 25}\\
a_7=\frac{-4}{\pi^2 49}
$$


How to filter noise from continuous path in image

So entropy filter sounded good, but it didn't really work because the actual signal is also considered high entropy, even after adjusting the window.

What I did was take a derivative of all of the data in the time axis (left to right). Then I capped the maximum value. Then I integrated the derivative.

deriv = data[i, 1:] - data[i, :-1]
deriv[deriv > max_] = max_

The derivitave is $ \frac{dy}{dx }$. In my mind dy is a vertical line segment, dx is a horizontal line segment. Putting them together makes a right triangle, where the function is the hypotenuse is f, the horizontal leg is dx and the vertical leg is dy. Thus $ tan^{-1} {\frac{dy}{dx}} $ is the angle of elevation for f. In my mind, putting a maximum on the derivative is like putting a maximum angle on f. It is like saying, "Cap the angle of elevation for f. If it goes above 150 degrees from the horizon, make it smaller."

It seems to work. It limits the maximum angle of elevation for the signal. You can see the original signal (in red) and the signal after derivative filter (in green).

Then I take original signal minus derivative filter applied to original signal. This leaves me with just the noise band, which I can amplify, blur, and then subtract out.


Baseline correction for continuous EEG signals

A baseline is generally an unwanted artifact that alter interpretation or computation on signals. The concept is not totally well-defined, and several different names coexist, with similar meaning:  drift, continuum, trend, background.

It is often described as a large-scale, slowly-varying reference level, that should be removed to compare, or adjust the signal with respect to a globally flat or  zero-reference. It can be seen as a generalized (non-constant) offset or bias  correction. When removed, it is generally simpler to further process the signal, quantify parameters, extract features.

In some domains, this smooth reference may also  include random noise.

Similar problems are answered here:

How to remove or filter the drift problem in measured Strain signal?
Help with automated baseline removal


How do I interpret an image of a spectrogram of an audio wind signal?

The almost steady colored part between 1 s and 35 s, and below about 30 Hz may be the sign of a   low-frequency, continuous noise that seems to be common to mild winds. There seems to be occasional gusts of wind around 12 s to 17 s or 22 s to 25 s. However, the needle aspect of the  spectrogram coefficients may suggest to modify the window, and the little information on frequencies above 256/512 Hz would suggest to look closer at the low frequencies, possibly with a linear frequency scale to better see details.


For eg. Using a length of 8 pilot samples can we detect a channel of 20 taps?

If your OFDM system is sensibly designed, then, no, you need a pilot on every subcarrier. They don't necessarily have to be sent in the same OFDM symbol.

"Sensible" OFDM design means that you choose your subcarrier spacing no smaller than you have to. The reason that you're doing OFDM in the first place is always that you have a frequency-selective channel and would like to divide that big problem into many smaller problems by dividing that frequency-selective channel into ideally flat subcarriers. You'd call the bandwidth in which you can reliably consider the channel flat coherence bandwidth. If you're narrower than that, everything is fine, if your subcarrier is wider, you run into frequency selective fading trouble.

Now, you can't choose an arbitrarily large number of subcarriers – remember, that number is directly proportional to the temporal length of your OFDM symbol, and if you make that longer than it needs to be, you risk it becoming longer than the coherence time. That means your channel changes during your usage of it, and you run into temporal fading trouble.

So, these are competing constraints: you want as many subcarriers as possible for flatness reasons, but you want as little subcarrier as possible for constantness reasons.

Therefore, a well-designed OFDM system will have as little correlation between the channel coefficients of neighboring channels as possible, because if neighboring channels were reliably similar, then you'd be using more subcarriers than necessary.

Therefore, that system will simply have as many different subcarriers as it has subcarriers. You hence need a pilot on every single one of these, at least once in time.

Luckily, temporal changes tend to be relatively benign if your symbols aren't too long: phases tend to linearly shift, so extrapolation in time direction tends to work OK. Thus, modern OFDM systems tend to have "all pilot" symbols at the beginning of "superframes", and then pilots sprinkled across the time/frequency plane over the remaining OFDM symbols in that frame. That's sufficient to interpolate both in time as well as in frequency domain between these positions¹.

Look at LTE, for example (source of that image):

The red dots are the "sprinkled" pilots, and the green and yellow ones are "pilot-only"/synchronization OFDM symbols. The blue ones actually carry data.

¹ warning: he who interpolates OFDM symbols linear shall be disappointed by the performance compared to what he might have gotten using sinc interpolation.


You can use forms of the Goertzel algorithm as a fixed length block filter, as a IIR filter with a very long impulse response, or as a fixed length running filter.

As a fixed length filter, a Goertzel results in a FIR filter similar (almost identical) to the result of 1 bin of a DFT or FFT of the same length.

As an IIR filter, a Goertzel filter is on the edge on instability (pole right on the unit circle), so you may need to bound the filter state and/or output externally, or reduce the gain slightly below 1.0.

As a fixed length running filter, you subtract the output of  2 offset (but stabilized) Goertzel IIR filters to get a FIR filter or 1 bin of a DFT (length equaling the offset) that updates every input sample, and with very few MACs per sample.


Obtain the frequency response of a time domain chirp signal

1) On the time axis I have wrong values: the time history lasts 45 s whereas in the spectrogram plot I see values from 0 to 2000

RTFM on surf.  You'll need to do something like surf(new_t,new_f,X) instead of just surf(X).  Your new_f should be something like

new_f = linspace(0,Fs/2,256);

NOTE: I don't have access to matlab, and am not sure of the orientation of X. You may have to swap the parameter locations of new_t and new_f on surf.

  4) I actually need to see only the frequency-amplitude plot, not the 3D representation, so how can I do that? At the moment I am using view([90 0]), but is there a way to plot directly the 2D plot?

Do something like:

plot(new_t, max(X));

You may have to do max(X') instead, depending on the orientation of X.

I'll need to think about answers to 2 & 3. Watch this space.


Determining Toy Helicopter Orientation From Image

For yaw it's sufficent to have three distinct data points fixed to the helicopter in a plane parallel to the rotor rotation: different colored leds for example. You could put them on three thin radial arms fixed to the landing gear. Maybe colored reflective materials would work, too, but you'd have to compensate for environmental light changes.

Then getting the yaw angle is simple.  Suppose the lights are red, green, and blue 120 degrees apart and ccw when viewed from the top.

Then in the image you will see these lights in various x-coordinate orders. As the helicopter rotates ccw fromt the position shown in the diagram:

G R B   // green and blue in front; red in back
R G B   // green in front; red and blue in back
R B G   // red and gree in front; blue in back
B R G   // etc...
B G R
G B R

(Of course at the borders of these regions two lights will merge to a single point.  These can be handled as special cases.)

So with 3 lights you have split all posible yaws into 6 x 60 degree segments.

You can further get down to 30 degree resolution by just comparing the left and right "gap" widths. If pointing directly away from the camera is azimuth theta zero degrees and the red light is on the nose (as in the picture), then the algoritm is:

Xr = -sin(theta)
Xg = -sin(theta + 60)
Xb = -sin(theta + 120)
if (Xg < Xr && Xr <= Xb) {
  // Green and blue are in front of red
  dxLeft = Xr - Xg
  dxRight = Xb - Xr
  yaw = (dxLeft < dxRight) ? 15+-15 : -15+-15;
}
else if (Xr < Xg && Xg < Xb) {
  // Green is in front of red and blue
  dxLeft = Xg - Xr
  dxRight = Xb - Xg
  yaw = (dxLeft < dxRight) ? 30+-15 : 60+-15;
}
else ... four more cases

The X+-15 above mean X degrees plus or minus 15: a 30 degree segment.  To get the exact offset, you can construct a preprocessed table that converts the ratio

R = dxLeft < dxRight ? xLeft / xRight : xRight / xLeft

into an an exact azimuth offset between -15 and 15.  This is pretty simple trigonometry to work out, and the nice thing is that it works independent of how far away the helicopter is, so long as it never turns upside down!

Another side benefit: the image distance between the outer two lights gives you inverse distance to the camera.


Your second plot is the square of the magnitude of the frequency response in dB, which indeed doesn't make much sense. You must use $|H(j\omega)|^2$, and not $(20\log_{10}|H(j\omega)|)^2$, as you did. Of course, $|H(j\omega)|^2$ will not give you any new information compared to $H(j\omega)$.


OK, I figured out what the problem was using a hex editor. The compiler on Windows was seemingly introducing random 0x0D throughout the file. That's because of the mode I was opening the file with, "w". fwrite converted every instance of '\n' into "\r\n". All I had to do was change the file opening mode to "wb+", and it sounds exactly like the linked audio file, after changing the carrier amplitude to 2e5.


Impulse response of a continuous system sampled with zero-order hold

People in ##dsp on Freenode helped me realize that the impuse response is a combination of impulse responses of two simple systems connected in series, that is, convolution of two impulse responses in the time domain.

$$
\begin{aligned}
G(z) &=
K \frac
  {\left( 1 - e^{-\frac{T_s}{T}} \right)}
  {\left( z - e^{-\frac{T_s}{T}} \right)}.
\\ &=
K \left( 1 - e^{-\frac{T_s}{T}} \right)
z^{-1}
\frac{1}{\left( 1 - e^{-\frac{T_s}{T}} z^{-1}\right)}
\end{aligned}
$$

$$
\begin{aligned}
h[n] &= \mathscr{Z}^{-1} \left\{ G(z) \right\}
\\ &=
K \left( 1 - e^{-\frac{T_s}{T}} \right)
\mathscr{Z}^{-1} \left\{
z^{-1}
\frac{1}{\left( 1 - e^{-\frac{T_s}{T}} z^{-1}\right)}
\right\}
\\ &=
K \left( 1 - e^{-\frac{T_s}{T}} \right)
\mathscr{Z}^{-1} \left\{
z^{-1}
\right\}
*
\mathscr{Z}^{-1} \left\{
\frac{1}{\left( 1 - e^{-\frac{T_s}{T}} z^{-1}\right)}
\right\}
\\ &=
K \left( 1 - e^{-\frac{T_s}{T}} \right)
\delta[n-1] * \left( e^{-\frac{T_s}{T}n} u[n] \right)
\\ &=
K \left( 1 - e^{-\frac{T_s}{T}} \right)
e^{-\frac{T_s}{T}(n-1)} u[n-1]
\end{aligned}
$$

The impulse response is exponential decay shifted to the right by 1.


Note the convention of notation: "$X[\cdot]$" is not the same as "$X(\cdot)$" (but they are defined solely in terms of the same "$x[\cdot]$").  The square brackets [] allow only discrete integer arguments, the round parenths () allow continuous real or complex arguments.

Discrete-Time Fourier Transform (DTFT):

$$ \mathscr{F}\big\{x[n]\big\} \triangleq X\big( e^{j\omega} \big) = \sum\limits_{n=-\infty}^{\infty} x[n] \, e^{-j\omega n} $$

This is literally the double-sided $\mathcal{Z}$-transform of $x[n]$ evaluated at $z=e^{j \omega}$.

If your original $x[n]$ has non-zero length of $N$, then DTFT is:

$$ X\big( e^{j\omega} \big) = \sum\limits_{n=0}^{N-1} x[n] \, e^{-j\omega n} $$

The length-$N$ Discrete Fourier Transform (DFT) of $x[n]$ is

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

The DTFT is a mathematical operation for the blackboard (or whiteboard or on paper).  The DFT is what computers do when they perform the Fast Fourier Transform, (FFT).

And the DFT is the DTFT of this finite-length sequence evaluated at $N$ equally-spaced frequencies between $0$ and $2\pi$.

$$ X[k] = X\big( e^{j\omega} \big) \bigg|_{\omega = 2 \pi k/N} $$

So now, what if you zero-pad $x[n]$ to a longer length $\tilde{N}>N$?

$$ \tilde{x}[n] = \begin{cases}
x[n] \qquad & 0 \le n < N \\
\\
0 \qquad & N \le n < \tilde{N} \\
\end{cases} $$

Then the DTFT of this is:

$$ \mathscr{F}\big\{\tilde{x}[n]\big\} = \tilde{X}\big( e^{j\omega} \big) = \sum\limits_{n=-\infty}^{\infty} \tilde{x}[n] \, e^{-j\omega n} $$

But the funny thing is that the two DTFTs are the same

$$ \tilde{X}\big( e^{j\omega} \big) = X\big( e^{j\omega} \big) $$

because there is no difference between $x[n]$ and $\tilde{x}[n]$ in the zero-extended DTFT.

But the DFTs come out a little different (only because $\tilde{N} \ne N$).

$$\begin{align}
\tilde{X}[k] &= \sum\limits_{n=0}^{\tilde{N}-1} \tilde{x}[n] \, e^{-j 2 \pi n k /\tilde{N}} \\
\\
&= \sum\limits_{n=0}^{N-1} \tilde{x}[n] \, e^{-j 2 \pi n k /\tilde{N}} \\
\end{align}$$

(Note the change of $\tilde{N}$ to $N$ in the second summation is because $\tilde{x}[n]=0$ for $N \le n < \tilde{N}$.  That has an interesting consequence below.)

Also this:

$$ \tilde{X}[k] = \tilde{X}\big( e^{j\omega} \big) \bigg|_{\omega = 2 \pi k/\tilde{N}} $$

Now it turns out that the inverse DFT is similar for $\tilde{x}[n]$

$$ \tilde{x}[n] = \tfrac1{\tilde{N}} \sum\limits_{k=0}^{\tilde{N}-1} \tilde{X}[k] \, e^{+j 2 \pi n k /\tilde{N}} $$

and the DTFT can be evaluated in terms of this expression for $\tilde{x}[n]$.

$$\begin{align}
 \tilde{X}\big( e^{j\omega} \big) &= \sum\limits_{n=0}^{\tilde{N}-1} \tilde{x}[n] \, e^{-j\omega n} \\
\\
&= \sum\limits_{n=0}^{N-1} \tilde{x}[n] \, e^{-j\omega n} \\
\\
 &= \sum\limits_{n=0}^{N-1} \tfrac1{\tilde{N}} \sum\limits_{ k=0}^{\tilde{N}-1} \tilde{X}[k] \, e^{+j 2 \pi n k /\tilde{N}} \, e^{-j\omega n} \\
\\
 &= \sum\limits_{k=0}^{\tilde{N}-1} \tilde{X}[k] \tfrac1{\tilde{N}} \sum\limits_{n=0}^{N-1}  e^{+j 2 \pi n k /\tilde{N}} \, e^{-j\omega n} \\
\end{align}$$

The DFT of $x[n]$ you seek is

$$\begin{align}
X[k] &= X\big( e^{j\omega} \big) \bigg|_{\omega = 2 \pi k/N} \\
\\
&= \tilde{X}\big( e^{j\omega} \big) \bigg|_{\omega = 2 \pi k/N} \\
\\
&= \sum\limits_{ \tilde{k}=0}^{\tilde{N}-1} \tilde{X}[\tilde{k}] \tfrac1{\tilde{N}} \sum\limits_{n=0}^{N-1}  e^{+j 2 \pi n \tilde{k} /\tilde{N}} \, e^{-j\omega n} \bigg|_{\omega = 2 \pi k/N} \\
\\
&= \sum\limits_{ \tilde{k}=0}^{\tilde{N}-1} \tilde{X}[\tilde{k}] \tfrac1{\tilde{N}} \sum\limits_{n=0}^{N-1}  e^{+j 2 \pi n \tilde{k} /\tilde{N}} \, e^{-j2 \pi n k/N} \\
\\
&= \sum\limits_{ \tilde{k}=0}^{\tilde{N}-1} \tilde{X}[\tilde{k}] \tfrac1{\tilde{N}} \sum\limits_{n=0}^{N-1}  e^{+j 2 \pi (\tilde{k}/\tilde{N} - k/N)n} \\
\end{align}$$

Note the $\tilde{k}$ is just a different "$k$" as a token counting variable in the summation.

There is this exponential summation identity in closed-form:

$$ \sum\limits_{m=0}^{M-1} a^m = \frac{a^M-1}{a-1}$$

so

$$\begin{align}
\tfrac1{\tilde{N}} \sum\limits_{n=0}^{N-1}  e^{+j 2 \pi (\tilde{k}/\tilde{N} - k/N)n}  &= \tfrac1{\tilde{N}} \sum\limits_{n=0}^{N-1}  \big(e^{+j 2 \pi (\tilde{k}/\tilde{N} - k/N)}\big)^n \\
\\
&= \frac{\big(e^{+j 2 \pi (\tilde{k}/\tilde{N} - k/N)}\big)^N - 1}{\tilde{N} \big(e^{+j 2 \pi (\tilde{k}/\tilde{N} - k/N)} - 1 \big)} \\
\\
&= \frac{e^{+j \pi (\tilde{k}N/\tilde{N} - k)}\big( e^{+j \pi (\tilde{k}N/\tilde{N} - k)} - e^{-j \pi (\tilde{k}N/\tilde{N} - k)} \big)}{\tilde{N} e^{+j \pi (\tilde{k}/\tilde{N} - k/N)}\big(e^{+j \pi (\tilde{k}/\tilde{N} - k/N)} - e^{-j \pi (\tilde{k}/\tilde{N} - k/N)} \big)} \\
\\
&= \frac{e^{+j \pi (\tilde{k}N/\tilde{N} - k)}e^{-j \pi (\tilde{k}/\tilde{N} - k/N)}\big( e^{+j \pi (\tilde{k}N/\tilde{N} - k)} - e^{-j \pi (\tilde{k}N/\tilde{N} - k)} \big)/(2j)}{\tilde{N} \big(e^{+j \pi (\tilde{k}/\tilde{N} - k/N)} - e^{-j \pi (\tilde{k}/\tilde{N} - k/N)} \big)/(2j)} \\
\\
&= \frac{e^{+j \pi (N-1)(\tilde{k}/\tilde{N} - k/N)} \sin\big( \pi N (\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)}{\tilde{N} \sin\big( \pi (\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)} \\
\end{align}$$

This shows the DFT of the original data (at the original size of $N$) can be extracted from the DFT of the zero-padded data:

$$\begin{align}
X[k] &= \sum\limits_{ \tilde{k}=0}^{\tilde{N}-1} \tilde{X}[\tilde{k}] \tfrac1{\tilde{N}} \sum\limits_{n=0}^{N-1}  e^{+j 2 \pi (\tilde{k}/\tilde{N} - k/N)n} \\
\\
&= \sum\limits_{ \tilde{k}=0}^{\tilde{N}-1} \tilde{X}[\tilde{k}] \left( e^{+j \pi (N-1)(\tilde{k}/\tilde{N} - k/N)} \frac{\sin\big( \pi N (\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)}{\tilde{N} \sin\big( \pi (\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)} \right)  \\
\end{align}$$

It ain't purty, but I think this is the simplest I can make it.  This is one form of expressing the Dirichlet kernel.

Edit note:

I made a completely legitimate change of $\tilde{N}$ to $N$ (because $\tilde{x}[n]=0$ for $N \le n < \tilde{N}$) and this results in some changes in the answer.  With $\tilde{N}$ unchanged the previous answer (which is just as correct) was:

$$ X[k] = \sum\limits_{ \tilde{k}=0}^{\tilde{N}-1} \tilde{X}[\tilde{k}] \left( e^{+j \pi (\tilde{N}-1)(\tilde{k}/\tilde{N} - k/N)}\frac{\sin\big( \pi \tilde{N}(\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)}{\tilde{N} \sin\big( \pi (\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)} \right) $$

They both have to be correct.  But clearly these two expressions are not equal:

$$ e^{+j \pi (\tilde{N}-1)(\tilde{k}/\tilde{N} - k/N)}\frac{\sin\big( \pi \tilde{N}(\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)}{\tilde{N} \sin\big( \pi (\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)}
\ne e^{+j \pi (N-1)(\tilde{k}/\tilde{N} - k/N)}\frac{\sin\big( \pi N (\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)}{\tilde{N} \sin\big( \pi (\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)} $$

So how is it that the summations (with the very same $\tilde{X}[\tilde{k}]$ coefficients) add up to be the same?

I would be very interested in anyone's answer to that.


Have a look at the examples for fir1, fir2, or firpm. The first two are straightforward, while firpm is a little more complicated as it derives the minimum filter length needed to match the design criteria, but it's handy to know.

For a simple lowpass filter this should be enough: Feel free to play around with N or look into using firpm to get the minimum order of N.

fc = 100; fs = 1024; N = 256;
b = fir1(N, 2*fc/fs);
freqz(b, 1, 65536, fs)


How to smoothen signal with missing values before differentiation?

The best tool for this job is normalized convolution. It can deal with missing samples as well as uncertainty.

The paper describing the method is "Normalized and Differential Convolution -- Methods for Interpolation and Filtering of Incomplete and Uncertain Data" by Hans Knutsson Carl-Fredrik Westin. There is a PDF on Semantic Scholar.

Normalized convolution basically takes a convolution with a Gaussian of the data, with 0 for the missing samples, and divides sample-wise with the result of the same convolution on a weight signal, which is 1 for known samples and 0 for the samples you set to zero. More formally:

$$
g = (( f m ) * h ) / (m * h)
$$

(With $h$ the convolution kernel and $m$ the weight signal). You can derive the analytic derivative of that equation, leading to a derivative operator that is robust to missing samples. See the paper for more details.


noises gets added therefore it is additive

correct!

which make me think that noises are not destructive

incorrect :(
simple thought experiment: You flip a fair coin $X$ (Head = -1 / Tail = 1) and tell me the result. The entropy here is 1 bit, i.e. the (expected) information ($I(X=\xi) = -\log_2 \left[P(X=\xi)\right]$) of each outcome is 1 bit.
Then there's additive noise $N$ that takes one of the values $\{-2,0,+2\}$ with equal probability.
When you receive a -1, you can't know whether the coin was Head and there was 0 noise, or the coin was Tail and there was -2 noise. Both are equally likely!¹
So, your additive noise is absolutely able to destroy information and hence is very destructive to your signal.
If you're more coming from a wireless communications background: Your $X\in\{-1,+1\}$ can be interpreted as BPSK. Now you see how even benign Gaussian noise destroys your reception when its sign is the opposite of your transmit symbol!

¹ we can even formalize that. Since $X$ (2 options) and $N$ (3 options) are independent, and each of them equidistributed, there's six possible combinations, each of them equally likely
 X | N  | Y = X+N
------------------
-1 | -2 | -3
-1 |  0 | -1
-1 | +2 | +1
+1 | -2 | -1
+1 |  0 | +1
+1 | +2 | +3

Thus, we have four possible outcomes for the sum of signal and additive noise, -3, -1, +1 and +3.

If we see +3 or -3, we get 1 bit of the 1 bit info in the coin toss. (It must have been +1, otherwise we couldn't get +3, or -1 for -3, respectively.) That happens in 2 out of 6 times, so with probability 1/3.
If we see -1, we don't know whether it was +1-2 or -1+2, so we have zero bit of the 1 bit of the coin toss. Same for +1. That happens 4 out of 6 times, so with probability of 2/3.

Thus, the expected information to get out of this channel is 1/3·1 + 2/3·0 bit = 1/3 bit, where put in full 1 bit! That's a very destructive additive noise channel.


The comment is mistaken.  The carrier frequency is 100 Hz.  The sample rate is 9600 Hz.  The carrier frequency divided by the sample rate is the carrier phase increment (phase change per sample), not the carrier frequency.  The units are radians / $2\pi$ (i.e. the phase increment is normalized from 0 to 1, not 0 to $2\pi$).  That phase increment can be used to index into a sine table to create a sinusoid.


Relationship between Bandwidth and Bit Rate

Assume you have a pulse of width $T$ which is your bit.

$$w(t)=1\quad for\quad-\frac T2\le t\le\frac T2\quad and\; 0\;elsewhere$$
 Your bit in frequency domain is:

$$W(f)=Sinc(T.f)$$

Look up Sinc function and you can see $Sinc(T.f)$ main lobe is between $\pm \frac 1T$. As your bit is widened in time, it occupies less bandwidth in frequency domain($T$ increases, $\frac 1T$ decreases) and vice versa. This is a very important result to understand and remember for the rest of your life. It's the reasoning behind relationship of frequency $BW$ and pulse duration $T$. Carrier frequency has nothing to do with it, it just shifts your entire signal spectrum up in frequency to make it suitable for transmission in the physical medium(channel) of your choice.

You should review Fourier transform and its properties.  Then look up Convolution theorem, which tells if you multiply signals in time, they get convolved in frequency domain. Then learn about Sinc function and its relationship to Rect function.


Mutual Information, is this calculation correct or should I contact the author about fixing it?

No there is nothing wrong with the formulas in the paper.

MI(S) is the average mutual information of set S over all classes c in C. It is NOT a reformulation of the formula for MI(W,C).
MI(w) is identical to your formula for I(Wi,C)  (except it uses slightly different notation).


Using fourier coefficients to reconstruct data in matlab

For real signals, it works a little differently than the complex case as there is a symmetry in the DFT that can be exploited.

The DFT and inverse DFT are conventionally defined as:

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

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

Now express your bin value in cartesian coordinates (due to the IMO poor convention of using x for the signal and X for the DFT bins, this is a different x):

$$ X[k] = x_k + i y_k $$

$$ X[N-k] = x_k - i y_k $$

Note:

$$ e^{i n (N - k) \frac{2\pi}{N} } = e^{-i n k \frac{2\pi}{N} }    $$

Evaluate the kth $(k>0)$ term:

$$
\begin{aligned}
T_k &= X[k] e^{i n k \frac{2\pi}{N} } + X[N-k] e^{-i n k \frac{2\pi}{N} } \\
&=  ( x_k + i y_k )( \cos( n k \frac{2\pi}{N} ) + i \sin( n k \frac{2\pi}{N} ) ) \\
&+  ( x_k - i y_k )( \cos( n k \frac{2\pi}{N} ) - i \sin( n k \frac{2\pi}{N} ) ) \\
&=  2 x_k \cos( n k \frac{2\pi}{N} ) - 2 y_k \sin( n k \frac{2\pi}{N} ) )
\end{aligned}
$$

Put that together, with a shift of one for MATLAB indexing fiasco, and you got what you want.

For the continuous reconstruction, then it is typical to rescale the domain from $ 0 \to N $ to $ 0 \to 2\pi $.

$$ t = \frac{n}{N} 2\pi $$

$$ n = t \frac{N}{2\pi} $$

The $k$th term then becomes:

$$ T_k = 2 x_k \cos( k t ) - 2 y_k \sin( k t ) $$

So for a real valued signal, the DFT real and imaginary are half the Fourier Series coefficients.  But you still have to divide them by N if you used the conventional unnormalized DFT.

Now, when we get to the Nyquist frequency when N is even, for real valued signals, to keep the results real:

$$ T_{N/2} = x_{N/2} \cos( \frac{N}{2} t ) $$

For odd N, the last $k$ will be (N-1)/2 and nothing special needs to be done.

Now when you add them up, you still need to divide by N, so in conclusion:

$$ a_k = \frac{2}{N} x_k $$

$$ b_k = -\frac{2}{N} y_k $$

Finally the zeroth bin, $k=0$, aka the DC bin, MATLAB index of one (stupid), for a real signal is the offset you start your series with.

$$ a_0  = \frac{1}{N} x_0 $$
$$ b_0 = 0 $$

For N even:

$$ a_{N/2}  = \frac{1}{N} x_{N/2} $$
$$ b_{N/2} = 0 $$

I know this isn't the clearest or cleanest explanation, but you can figure it out.

Ced


What is Earth Movers Distance and how is it used in pattern recognition?

A set of clusters. Each cluster is represented by a pair with (1) one single point (usually the center in the cluster) ,and (2) a weight coefficient that usually in proportional to the size of the corresponding cluster. Note that the point can be anything including the coordinates of the pixel, the color information, intensity...etc.


Running a sox spectrogram on a bandpass filtered wav file

sox is doing exactly what you tell it to do – it plots a spectrogram of the bandwidth you can represent with your sampling rate. The fact that most of your picture shows no energy just shows your filter is working, and that's exactly what one would want to see!

However, as you can see, your signal is solidly oversampled: You can downsample it to a much lower sampling rate after filtering! (and downsampling by a factor $N$ here really is just throwing away $N-1$ samples, keeping one, throwing away $N-1$ samples..., since your filter is an excellent anti-aliasing filter.)


Filter Design near Nyquist frequency (complex curve fitting)

You can't get a good fit very near the Nyquist frequency (or near DC) because information is lost due to the imaginary components around that frequency being of (almost) opposite phase and thus (nearly) cancelling out.  Thus, there will be ambiguity errors due to trying to fit any curve to lost information.

If possible, you can try to gather more information near (but not at) the Nyquist frequency by sampling for a longer time span or at a higher sample rate.  This will move the measurement point slightly farther from where there is complete cancellation of the imaginary component.


How to prove that frequency mixer is a nonlinear system?

If the system is nonlinear then if $y_1(t)$ is the response to the signal $x_1(t)$, and  $y_2(t)$ is the output given input signal $x_2(t)$ then the response to the signal

$$x(t)=a_1x_1(t)+a_2x_2(t)\tag{1}$$

with arbitrary constants $a_1$ and $a_2$ will generally not be equal to

$$y(t)=a_1y_1(t)+a_2y_2(t)\tag{2}$$

However, for the given system an input signal (1) always results in an output signal (2). The conclusion is that the system is linear.

The confusion arises because you think that linear systems cannot add new frequencies to the input signal. This is generally not true. Only linear and time-invariant (LTI) systems have that property. Linear but time-varying systems can produce new frequencies in the output signal, and a modulator is a classic example of a linear time-varying system.


It does the same thing (i.e. shift the frequency to another band), though mathematically it differs a bit. Fourier transforms of the two choices are as shown below:

$$
\begin{array}{|c|c|}\hline
f(x)\cos(ax) & \displaystyle\frac{\hat{f}\left(\xi-\frac{a}{2\pi}\right)+\hat{f}\left(\xi+\frac{a}{2\pi}\right)}{2}\\
\hline
f(x)\sin(ax) & \displaystyle\frac{\hat{f}\left(\xi-\frac{a}{2\pi}\right)-\hat{f}\left(\xi+\frac{a}{2\pi}\right)}{2i}\\\hline
\end{array}
$$

The negative frequencies not only get shifted, but also flipped in polarity

Now for the derivation of these formulas:
Using the Great Mr. Euler's identity: $e^{iax} = \cos(ax) + i \sin(ax)$

$$\cos(ax) = \frac{e^{iax} + e^{-iax}}{2}$$ and
$$\sin(ax) = \frac{e^{iax} - e^{-iax}}{2i}$$

\begin{align}
\textrm{F.T.}\left\{f(x) \cos(ax)\right\} &= \int\limits_{-\infty}^{\infty} f(x) \cos(ax) e^{-i2\pi\xi x} dx\\
&= \int\limits_{-\infty}^{\infty} f(x) \frac{e^{iax} + e^{-iax}}{2} e^{-i2\pi\xi x} dx\\
&= \frac{\int\limits_{-\infty}^{\infty} f(x) e^{-i2\pi(\xi + \frac{a}{2\pi} )x} dx + \int\limits_{-\infty}^{\infty} f(x) e^{-i2\pi(\xi - \frac{a}{2\pi})x} dx}{2}\\
&= \frac{\hat{f}(\xi + \frac{a}{2\pi} ) + \hat{f}(\xi - \frac{a}{2\pi} )}{2}
\end{align}

You can derive the formula for $\sin$.