Quadrature signal acquisition simulation

For a causal time domain signal, the Fourier Transform will be complex with the imaginary part as the Hilbert Transform of the real part.
This may be clearer by noting the following additional properties of the Fourier Transform:

The FT of a non-causal real even waveform ($f(t) = f(-t)$) will be
all real.
The FT of a non-causal real odd waveform ($f(t) =
   -f(-t)$) will be all imaginary.
The FT of a non-causal imaginary odd waveforrm will be all real.
The FT of a non-causal imaginary even waveform will be all imaginary.

It is easy to prove all of the above by observing waveforms as rotating phasors of constant magnitude in time on a complex IQ plane, since the FT is the mangitude of each of those phasors at the position in frequency based on their rate of rotation. (In fact the Fourier Series Expansion is specifically to decompose any arbitrary single valued analytic function into such rotating phasors).
With those properties in mind, you can extend that to causal and anti-causal functions by adding or subtracting the cases listed above: Causal and anti-causal functions are the sum of an even and odd functions. Understanding this provides more insight into many other Fourier Transform properties.
As far as the first statement given about the interesting relationship between the real and imaginary components; if it is still not clear how this occurs it may be easier for some to swap the time and frequency domains and consider what occurs: In the time domain if you add the Hilbert Transform of a real signal as an imaginary component the result is a one-sided spectrum. A one-sided spectrum is the equivalent in the frequency domain to what a causal or anti-causal signal is in the time domain. Stated another way, a causal signal is a one-sided time-domain signal.
Consider the very simple case of a cosine which is a two sided spectrum with a positive negative frequency. This is also apparent in Euler's Identity where we see the two rotating phasors I was mentioning above, one rotating counter-clockwise (positive frequency) and one rotating clockwise (negative frequency):
$$cos(\omega t) = \frac{1}{2}e^{j\omega t} + \frac{1}{2}e^{-j\omega t}$$
Now consider the identity given as:
$$e^{j\omega t} = cos(\omega t) + j sin(\omega t)$$
$sin(\omega t)$ is the Hilbert Transform of $cos(\omega t)$, and $e^{j\omega t}$ is the single sided spectrum referred to earlier.
So we see that the following relationship holds:
$$f(t) = a(t) + j H\{a(t)\}$$
Where $f(t)$ is a function with a single-sided Fourier Transform, and $H\{\}$ is the Hilbert Transform.
Similarly by swapping time and frequency domain the following relationship will hold:
$$F(\omega) = A(\omega) + j H\{A(\omega)\}$$
Where here $F(\omega)$ is the Fourier Transform of a causal time domain function $f(t)$.  The resulting $f(t)$ which has a single-sided Fourier Transform is referred to as the "analytic signal" and it also has the property that it will be minimum phase: it's Laplace transform will have all poles and zeros in the left half plane (and similarly for discrete-time systems the z-transform will have all poles and zeros inside the unit circle). A minimum phase system will have a single-sided frequency spectrum but not all single-sided frequency spectrums are minimum phase systems: We can cascade a minimum phase system with an all-pass filter that only modifies phase resulting in a mixed-phase or maximum-phase system, but since this only modifies phase in the spectrum and doesn't add new frequencies, the spectrum will still be single-sided.
For further details on this see:
FFTs of a complex signal - separating the real and imaginary parts


Mathematical model of a signal in Compressive Sampling

The model is the Fourier series of a signal up to the W/2th component, with sparse coefficients. That is, most of the set $\{a_k\}_{k=-W/2}^{W/2}$ is zero.

It is because the signal is bandlimited, and thus for a real signal, $a_{W/2}$ will be real. However, if this constraint is not met it is better to use $\pm W/2$
For a real valued signal, the coefficient of the negative frequency is just the conjugate of the positive frequency coefficient. It really just is a mathematical artifact. You could also write the model as

$$f(t) = \sum_{\omega\in\Omega} b_w \cos(2\pi \omega + \phi_\omega)$$

where $b_\omega$ is real and $\phi_\omega \in [0,2\pi)$ is a real phase value.

That is, the model represents the signal as the sum of sinusoids of different amplitudes and phases.

It will be the Fourier series approximation of the signal using just the first two frequencies.

You should take a look at this paper and the example code:

Towards a mathematical theory of super-resolution


If underwater communications is as complicated as I expect it to be, you get reflections from your own vessel, the ground, larger swarms of fish, and doppler due to water currents, and interesting pressure-dependent nonlinearties. An interesting field, indeed!
Now, that puts your modem in the position to have to compensate a lot of interesting¹ effects:

High attenuation
Doppler due to currents
Reflection on you vessel, ground, fish, submarines, cables, density strata, and hence a multipath channel
Frequency drift
potentially strong nonlinearity due to transducers and density effects in greater depths
…
Temporal variation of all the above effects

Thus, your modem will have to not only add error correction to the data bits you give to it, map these bits to some sort of symbol, and modulate that, and reverse that at the receiver. It will also have to estimate the channel and revert it; that's what frequency, phase and timing synchronization are all about, and what equalization does to revert a multipath channel as far as possible.
But the channel estimation would fall flat the moment multiple modems are transmitting.
Hence, you can't just add CDMA as extra layer to this problem that you can address on the data-providing side of the modem. You need to redesign the modem to be able to work as CDMA modem.

How to Tx CDMA signal physically using an existing point-to-point modem?

Redesign the modem.
Since you say in the comments that you've worked with SDRs: exactly like that, you'd develop your whole system in SDR software, simulate it, evaluate performance and improve it before you start working with "real" hardware transducers. One of such might be your current modem, but with a pretty different firmware that would be a port of your prototype to the specific Blackfin platform.
Maybe this talk (slides, paper) is a nice starting point for you; they do long reach, low speed communications using SDR approaches.

¹ and "interesting" here is the engineering definition of the word: "your rocket explodes whenever someone installs a printer in the office building? *Interesting!*"


Are ACF plots sufficient to detect white noise assuming probable existence of non linear associations?

In other words, are ACF plots sufficient for determining whether a signal is a white noise?

A plot is never sufficient. An impulse-shaped ACF is. (A plot can only ever be an estimate based on observation; to find the true ACF, you'd have to observe the process for infinity, which you can't.) You can only observe that the plot approximates the ideal ACF sufficiently much (to your definition of "sufficient").
Anyway, assuming you mean "ACF being a dirac delta", yeah, that fulfills the definition of white. That is sufficient, because it literally is the definition.

can there be a scenario where an ACF plot that is considered to reflect white noise be actually representative of non-linear associations

No, because white noise is defined to be any random process with a delta dirac ACF.

certain non-linear associations can exhibit zero autocorrelations

Yep! "associations" is a bit of very vague term, though. I'd prefer if you said:

Are there discrete random processes (i.e. random sequences) where elements aren't independent from each other, but still be uncorrelated?

Independence is a strictly defined term: It means that the joint density functions of multiple points in time is the product of the individual densities at these times. I.e.: Let $X(t, \xi)\in \mathbb C, \, t \in \mathbb Z$ be a process at times $t$ (with different realisations $\xi$).
Then, the ACF is defined as
$$\phi_{XX}(t_1,t_2) = \mathbb E\left(
X(t_1)X^*(t_2)
\right)$$
and a process is white iff (assuming wide-sense stationarity, so that $\phi(t, t+\tau) =: \phi(\tau) \forall t$)
$$\phi_{XX}(t_1,t_2)=\phi_{XX}(\tau)=\begin{cases}
> 0 & t_1=t_2\iff \tau=0\\
= 0 & t_1\ne t_2\iff \tau \ne 0
\end{cases},$$
from which it's easy to see that different elements from an observed series are uncorrelated.
However, as you probably know, correlation is not the same as independence.
Counterexample:
$$Y(t) = \begin{cases}
U(t) & t\ge 0\\
S\cdot Y(-t) & t < 0
\end{cases},$$
with $U(t)\sim\mathcal U([-1,1])$ (independently uniformly drawn from the $[-1,1]$ interval), and $S\sim\mathcal U(\{-1,1\})$ (a randomly independently drawn sign).
Obviously, different elements from the sequence are not independent. $Y(-5)$ is either $Y(5)$ or $-Y(5)$. Other values are not allowed.
Checking the ACF, you'll see
\begin{align}
\phi_{XX}(t_1,t_2) &= \mathbb E\left(
X(t_1)X^*(t_2)
\right)\\
&=\begin{cases}
\mathbb E\left(
U^2
\right) =\frac13&\text{for }t_1=t_2&\text{(Variance of uniform continous RV)}\\
0 &\text{for } t_1\ne -t_2 &\text{(the $U(t)$ are uncorr.)}\\
\mathbb E\left(
S\cdot U^2
\right)
 &\text{for } -t_1=t_2
\end{cases}\\
&=\begin{cases}
\frac13&\text{for }t_1=t_2\\
0 &\text{for } t_1\ne -t_2\\
\mathbb E(S)\mathbb E\left(
U^2
\right)
 &\text{for } -t_1=t_2 & \text{(Since $S(t)$,$U(t)$ indep.)}
\end{cases}\\
&=\begin{cases}
\frac13&\text{for }t_1=t_2\\
0 &\text{for } t_1\ne -t_2\\
0\cdot\frac13
 &\text{for } -t_1=t_2 \\
\end{cases}\\
&=\begin{cases}
\frac13&\text{for }t_1=t_2\\
0 &\text{for } t_1\ne t_2
\end{cases},
\end{align}
which seems to say "this process is white".
BUT: this process is not Wide-Sense Stationary, so applying the term "white" to it doesn't "work".  So, you need to be very careful to what you apply the label "white", to begin with.


Restore RGB Values from Luminance Value?

No, because luminance is a weighted sum of RGB components and has no color information, so it is greyscale only.


calculating estimate of a state for a system with two observations of the state from different times

HINT

$$ \left| \begin{array}{cc} A & 0 \\ 0 & 1 \end{array} \right| $$


The sampling is indeed analogous to mixing as to my understanding. In the sampling process, we multiply the time domain signal with an impulse train  - the impulses in time are represented as impulses in frequency at integer multiples of the sampling rate. So instead of one or two (for a real sine wave) impulses in frequency, we have an infinite number but the process of multiplying in time is otherwise identical. The result you would get (which explains the aliasing well) is the same thing you would get if you had an infinite number of mixers and LO's, one for each harmonic of the sampling clock.

I have explained this with the additional graphics below to help those more familiar with RF mixing to understand sampling and aliasing, and then undersampling as well.

The reference to third order intermodulation product is not what would explain this aliasing, as that is specifically due to non-linearities in the signal chain. This can certainly occur for the same reason in ADC's causing intermodulation distortion but that is not what causes the aliasing.  So when you see other spurious products in your spectrum, these certainly can be due to non-linear distortion (which can be confirmed by modifying the power level of your dominant signal to see if these products change), or can be caused by aliasing from other frequency bands due to insufficient front end filtering prior to ADC conversion (which means they would be present regardless of the presence of your input signal or not). Another source is spurs on the sampling clock itself, which too is explained well as a mixing process.

To understand the aliasing mechanism and how it is identical to mixing, first observe the sampling process for a 3 Hz Cosine Wave as demonstrated in the graphic below. Note that since we multiply in time the cosine wave with the time domain impulses, we convolve in frequency the 2 impulses representing the cosine wave as in the top portion of the graphic below with the impulses in frequency represented by teh middle portion of the graphic, resulting in the digital spectrum as given in the bottom graphic. In this case the signal was sampled at 20 Hz, so the output spectrum repeats every 20 Hz, so really only the spectrum from -10 Hz to +10 Hz needs to be given for the digital spectrum, as shaded in the graphic. (Or from 0 to 20 Hz, basically any 20 Hz will completely represent the output spectrum). However when working with mixed signal or multi-rate systems, I find that it often mentally helps to "un-roll" the digital spectrum and represent it out to $+/- \infty$ as I have done in this graphic.

Note too that the output spectrum can be completely explained as a mixing process: The two impulses in the middle portion of the graphic at +20 Hz and -20 Hz do represent a real sinusoid at 20 Hz, while the two at +/- 40 Hz represent a sinusoid at 40 Hz, etc... Each of the outputs shown in the digital spectrum can be explained using the traditional "sum and difference" frequency outputs that you may be accustomed to when working with an RF mixer (as explained by the multiplication of two real sinusoids).

This explains aliasing quite well, as introduced in the undersampling example shown in the graphic below (Even though I explained the above with real signals and sinusoids, I prefer to work with complex frequencies such as $e^{j\omega t}$ recognizing that a cosine is $e^({j\omega t} + e^{-j\omega t})/2$ and from that we see that each impulse in these frequency domain plots is a single instance of $e^({j\omega t}$. Then instead of dealing with sums and differences, we just add the frequency terms. For example, the graphic below shows that the frequency component that is labeled "IF" results from the time domain multiplication of the input frequency component at "RF" with the specific frequency impulse in the impulse train that is labeled "-LO".

With that view we can easily see how aliasing occurs (as explained by "mixing" operations!), and the next two graphics show the identical result between sampling in the 1st Nyquist Zone vs 3rd Nyquist zone with the first graphic colored to show which frequency components of the digital sampling spectrum are responsible for which outputs in the digital output spectrum.


First, shift: put the minimum to $0$, by compensating the actual minimum $m=−18.3667⋅10^5$ for every pixel: $p\to p - m$. Now your pixels are between $0$ and a new maximum $M = 9.3127 - m$. Finally you want the final image in $[0,255]$. The second operation is scale: multiply by something, so that $0$ remains at $0$, and $M$ is cast to $255$. So you have to multiply by $255/M$.


Filtering $ n \times n $ Images by Separable $ m \times m $ Filters. Computation Time for Filtering Using FFT, 2D Convolution and Two 1D Convolutions

I'm not sure exactly what you're after, but just to try to add data:

Using the separability property of the filter is always the right choice.
Given it is separable we now have to apply 1D convolutions. Selecting which method to do it (Frequency Domain / Spatial Domain) depends on the length of the filter and the number of pixels in each dimension of the image.
There are cases where it is better to do FFT on the Rows and Spatial Convolution on the Columns.
This is all true when "Counting" FLOPS. Real life timing is more than that. for instance, if you're using highly tuned Convolution implementation and yet "Classic" DFT implementation you might be faster doing the Spatial way even for dimensions the analytic calculation predict it is better to be done in frequency domain.


As the commenters have stated, this is not the same. The relationships between down/up-sampling and filtering are formally described by the noble identities.

In your case down-sampling followed by filtering with [-1 0 1] is equivalent to filtering with [-1 0 0 0 1] first and than downsampling.

In any case you need to make sure your don't alias before down-sampling the picture.


Two voice pronunciation comparison similarity MFCC + DTW

I will answer your questions in reverse order.

4: DTW (Dynamic Time Warping) is not a library but an algorithm. It allows aligning two sequences by warping them in time. You can use it for pretty much any kind of sequential data, for which a metric (distance) can be defined. Generally, you calculate a distance between each and every point of both sequences and produce a big matrix of accumulated distances. In the end, a path that produces the smallest final cost is selected as the optimal.

Here is an example of two 1D signals: $\color{blue} {x_1}$ is an exponential sweep sine and $\color{orange} {x_2}$ is a sinusoid with a constant frequency. The goal is to warp the $\color{orange} {x_2}$ in time so it would match the $\color{blue} {x_1}$ as close as possible. The figure below shows the outcome of such alignment. The signal $\color{orange} {x_2}$ had been stretched so that it matches $\color{blue} {x_1}$. However, $\color{orange} {x_2}$ is shorter than $\color{blue} {x_1}$, hence it had to be padded towards the end. The subplot at the bottom shows the evolution of the distance for each frame, where the final value is the accumulated cost - the very thing that you are using in your comparison. From this example, you can immediately see some limitations of the DTW.

Just to wrap up, here is the full DTW matrix together with the optimal path.

3 & 2: It's a distance, so the lower the better. The lowest possible distance is $0$, which can be achieved when:

Both signals are exactly the same
One signal is the stretched version of the other one, such as [0, -1, 0, 1, 0] and [0, -1, 0, 0, 1, 1, 0].

1: As I've mentioned before, DTW compares the distance between each and every point and then accumulates the cost of all operations that lead up to that point - the path that will yield the minimal cost corresponds to the best alignment.

In terms of similarity, signals will be considered close to each other if stretching or squeezing certain parts will minimize the difference between the MFCC vectors for each frame.

These answers might be relevant for you: click, click, click.


Solving LASSO ($ {L}_{1} $ Regularized Least Squares) with Gradient Descent

It can easily solved by the Gradient Descent Framework with one adjustment in order to take care of the $ {L}_{1} $ norm term.

Since the $ {L}_{1} $ norm isn't smooth you need to use the concept of Sub Gradient / Sub Derivative.
When you integrate Sub Gradient instead of Gradient into the Gradient Descent Method it becomes the Sub Gradient Method.

In the case of the $ {L}_{1} $ norm the Sub Gradient is given by the function $ \operatorname{sign} \left( \cdot \right) $.
In this context it works on vector in Element Wise manner.

So the Pseudo Code (MATLAB) would look something like:

for ii = 1:numIterations
    vG = (mA.' * (mA * vX - vY)) + (paramLambda * sign(vX)); %<! Sub Gradient
    vX = vX - (stepSize * vG); %<! Sub Gradient Update
end

Yet, As mentioned by @David, this is not a modern way to solve this problem.
Iteration wise it seems that Coordinate Descent is the fastest method for solving the LASSO ($ {L}_{1} $ Regularized Least Square) problem.

You may look at a project I created (${L}_{1}$ Regularized Least Squares - Solvers Analysis) which compares many method for that optimization problem:

Implementation is included so you can see how it works.


If you look around for "image patch self-similarity" there are some things that are along those lines, for example the work of Eli Shechtman and Michal Irani

http://www.wisdom.weizmann.ac.il/~vision/SelfSimilarities.html

or the Calvin research centre:

http://www.vision.ee.ethz.ch/~calvin/software.html

(see Global and Efficient Self-Similarity for Object Classification and Detection)


I'm assuming this is the uplink channel since you want the received power to be different for each user.

A DSSS has a chip sequence that modulates the incoming signal. Basically the chips are just +1 and -1 and hence your TX signal power still is determined only by the incoming signal.

Generate random bits - A simple function randint.
Modulate them however you want - You can use qammod. - Check out the book "Fundamentals of Digital Communications" by Proakis - Chapter 3 for intricacies regarding energy Normalization for constellations.
Pass them through a Pulse Shaping filter of sufficient length. Check this link out for the same: http://www.mathworks.com/help/dsp/examples/pulse-shaping-filter-design.html
Multiply with the Chip Sequence. The way DSSS works is each user will have a Spreading code that is almost perfectly orthogonal to each other user. I think MATLAB can generate spreading sequences you might want to see: http://www.mathworks.com/help/comm/examples/spreading-sequences.html
Once you've generated the signals for all users, pass them each through a Flat fading Channel(you can consider different average received energy and generate the fade vector with variances appropriately - this will give you different SINRs) - You'd want your Chip period Tc to be very small compared to the DelaySpread of the channel and you can assume this which will hold the Flat Fading condition to be true. Multiply the fade vector for each user symbol by symbol with their corresponding signal points.
Receiver receives the Linearly combined signal and then you need to generate an equivalent length AWGN vector considering the average TX energy as the average TX energy you generated at the each user but the signals from each user suffered independent and different fades with different paths leading to different SINRs.
or AN ALTERNATIVE FOR (5,6): After spreading the sequences, alternatively you can multiply the whole signal for each user by a constant at the Transmitters before the fading channel leading to different energies(equivalent to different Path Losses for each signal) and linearly combine them all, Transmit and Receive on just one Fading Channel. Further you'd need to add AWGN considering the Average Signal energy of the signal you generated at the output of the Modulator(essentially the average energy of the Constellation).
The receiver now is like a Multi-User detector; you need to create copies of the received signal(as many number of users) and correlate with the Spreading sequences for each user.  You will essentially now remove the spreading code(Since synchronization is perfect) and you'll have the received signals for each user with different SINRs and so you can proceed to Matched Filtering using the same Filter you used for Pulse shaping at the Transmitter, Sampling, Demodulation using qamdemod and Detection using Maximal Likelihood(Minimum Distance).

You can refer to "Wireless Communications" by Andreas Goldsmith on further technical details in Chapter 3 for Fading and Chapter 13 and 15 for Multiuser Spread Spectrum Systems.

Note: Since you've titled this as Imperfect Power Control, I'd say the alternatives for 5,6 make more senes - the TXs are having bad knowledge about the channel and are not transmitting according to proper Powers.


Cosine of amplitude 1 but GNU Radio shows FFT has amplitude > 30,000

Consider the formula for the DFT (which the FFT efficiently computes as an algorithm):
$$X(k) = \sum_{n=0}^{N-1}x(n)e^{-j2\pi nk/N}$$
Notice that it is a summation over $N$ samples total. Also note using Euler's formula that a cosine function can be expressed as two exponential terms (that are visible in the DFT result) as:
$$x[n] = cos(2\pi f n + \theta) = \frac{1}{2}e^{j(2\pi f n + \theta)} + \frac{1}{2}e^{-j(2\pi f n + \theta)}$$
Observe how in the DFT equation that we multiply the two exponential terms in x[n] by $e^{-j2\pi nk/N}$, that latter which is a phase rotation (the exponents which is the phase add in the product, and generally the term $e^{j\phi}$ is simply a complex vector with magnitude $1$ and angle $\phi$ just in case there was any confusion in that; $e^{j\phi}$ is simply a phase rotation). For each of those two exponential terms, there is a value of $k$ that would cause one of them to rotate to the same phase for every $n$, which then would grow to the large value in the summation. Given the OP's summation has grown to approximately 34,000, and that the magnitude of each of the exponential terms is $1/2$, suggests that there are approximately 68,000 total samples in the DFT assuming no windowing is also applied to the signal (which would then require even more samples due to the amplitude reduction through windowing). For every bin in the DFT result, this is the correlation to $e^{j\omega_k n}$ (a spinning phasor at frequency $\omega_k$, with $\omega_k = 2\pi k/N$), and sinusoidal functions have two such terms as given my Euler's formula: The worlds should now align.
$f$ is the fractional frequency term typically given as a value between $0$ and $0.5$ for the form in the equation above, where it represents the continuous time frequency divided by the sampling rate. Here we see further the values for k where that occurs:
When $k = Nf$
$$X(k=Nf) =  \sum_{n=0}^{N-1}x(n)e^{-j2\pi fn}$$
$$ = \sum_{n=0}^{N-1}\bigg( \frac{1}{2}e^{j(2\pi f n + \theta)} + \frac{1}{2}e^{-j(2\pi f n + \theta)} \bigg)e^{-j2\pi fn}$$
$$ = \frac{1}{2}\sum_{n=0}^{N-1}e^{j\theta} + \frac{1}{2}\sum_{n=0}^{N-1}e^{-j(4\pi fn + \theta)}$$
The first term grows to $\frac{N}{2}e^{j\theta}$ in the summation while most of the individual phasors in the second term will cancel leaving either zero or small residual depending on the value for $f$.
Similarly consider when $k=N(1-f)$
$$X(k=N(1-f)) =  \sum_{n=0}^{N-1}x(n)e^{-j2\pi(1-f)n}$$
$$ = \sum_{n=0}^{N-1}x(n)e^{-j2\pi n}e^{+j2\pi fn} = \sum_{n=0}^{N-1}x(n)e^{+j2\pi fn}$$
$$ = \sum_{n=0}^{N-1}\bigg( \frac{1}{2}e^{j(2\pi f n + \theta)} + \frac{1}{2}e^{-j(2\pi f n + \theta)} \bigg)e^{+j2\pi fn}  $$
$$ = \frac{1}{2}\sum_{n=0}^{N-1}e^{4\pi fn + \theta} + \frac{1}{2}\sum_{n=0}^{N-1}e^{-j\theta}$$


Can we have a Digital Anti Aliasing filter?

I agree with pichenettes's answer but I would like to add that it is pretty common practice to use a simple inexpensive low-order analog anti-aliasing filter, and do the rest of the anti-aliasing filtering in the digital domain. This implies of course that you do not process at maximum sampling rate but that you downsample after the digital anti-aliasing filter. Summarizing:

Of course you need an analog anti-aliasing filter.
The analog filter can be kept very simple if you can downsample your signal. In this case you can do more aliasing removal in the digital domain (before downsampling).


Side lobes not visible in FFT of rect function

You have to imagine your signal as being periodically extended forever in both directions. Since you square window is constant over the whole interval when it is periodically extended it becomes flat everywhere. In other words, you don't have a square window. You just have a flat line. That's why only the zero frequency band on the fft has any amplitude. Try doing it again with some zeros on either side of your square window.


AM modulated signal not centered at carrier frequency

Here is a code snipped that works, you can get it and modify according to your needs.

Fsim = 5E6;         % simulation sampling frequency
Tsim = 1/Fsim;      % simulation sampling period
Tdur = 0.1;         % duration of simulation in seconds...
t = [0:Tsim:Tdur];  % t = n*Tsim; sampled time instants...
L = length(t);      % number of samples in the signal

f_car = 1E6;        % DSB-SC modulation carrier frequency

bl_msg_time = filter(fir1(264,0.03),1,randn(1,L));   % generate a bandlimited message signal

msg_dsb_time = bl_msg_time .* cos(2*pi*f_car*t);
msg_dsb_freq = fftshift(fft(msg_dsb_time));
msg_dsb_freq=msg_dsb_freq/length(msg_dsb_freq);

fx = linspace(-Fsim/2,Fsim/2,L);
plot(fx,abs(msg_dsb_freq));
title(['DSB-SC modulated signal, Fcar = ' , num2str(f_car)])
xlabel('frequency in Hz')

The produced output is:

As you can see the modulated signal is centered at the correct frequency.


Power Spectral Density (PSD) from autocorrelation

I am not sure I follow your exact questions.

  Is there a specific mathematical reason as to why the latter method
  appears to carry more energy at the higher frequencies?

No. The two methods are equivalent and will match to the numerical precision you are using.

  Auto-correlation function is analogous to a low-pass filter, meaning
  low frequencies are getting through, and high frequencies are
  attenuated out (gradual loss) at a certain cutoff frequency.

I do not understand this claim. How is auto-correlation a low-pass filter? Auto-Correlation is auto-correlation.

  Technically both methods shows some form of aliasing. So are there
  different kinds of aliasing? One is observed to aliased over a range
  of wavenumbers, while the other shows a spike.

Your input signal is already sampled and you are not changing the sampling rate, so I do not understand this either.

If the PSDs as estimated from two methods do not match, then there is a mistake in your code. Below is an example using octave/matlab code:

N = 1000;
x = rand(N,1);
Nconv = 2*N - 1;

% PSD from autocorrelation.
a = conv(x,flipud(conj(x)));
A = abs( fftshift(fft(a)) );

% PSD from FFT.
X = fftshift(fft(x,Nconv));
A2 = X .* conj(X);

% The residual between the two estimates.
err = 10*log10( norm(A - A2) / norm(A) );


Removing a darkened caption with text on it in an image

Does this algorithm even make sense.

Your caption box is actually pretty wide – the algorithm that you've linked to would need to replace the complete content by its neighboring pixels. That won't look right – it will look as if someone "stretched" the guitar there.

So if you ask me, no, this is not the right approach here.

  Now what I noticed is that the caption is not fully opaque.

I was confused for a while until I realized that with "caption" you mean the box, and I mean the text. So I'll refer to the box to avoid further confusion.

So, yes, that would be what I'd look at here. The image just seems to be darkened in the box area – what about you just take this part of the image, and reverse that darkening, by increasing the brightness? I'd assume that degree of darkening is static – so a bit of trying out in e.g. GIMP you should be able to figure this out; I tried this, and  got this:

as you can see, it's slightly blueish- although the brightness seems to match. That might either be because whatever draws the box isn't actually reducing the brightness only, but somehow also changes the hue of the image in the region, but it might very well also be an artifact of JPEG compression; notice that your picture has a terrible quality to begin with. Or this is a color space conversion effect.


how to visualize image segmentation using matlab

some pointers for you to figure out the solution:

you can convert segmentation results to CONTOURS and then plot contours, with thick boundary lines.
first convert the image to grayscale,

x = rgb2gray(x);
img(:,:,1) = x;
img(:,:,2) = x;
img(:,:,3) = x;

now, for each cluster you can modify the color for that part of image. say for example row1, col1, keep the index for class 1:
  for i=1:size(row1,1)
        img2(row1(i), col1(i), 1) = 200;
    end
in the end show image with modified colors.
show image, then hold on, then show ONLY one cluster at a time as a region filled with a solid color, then set transparency to 50%.

two ways

to set transparencey:
set(h, 'AlphaData', alpha_data); % e.g. alpha_data = 0.5
or
h=imagesc(img2);
alpha(h,0.5);


Why are there two phases for one dominant frequency?

If the dominant frequency is not exactly integer periodic in the FFT width, then most of that frequency peak's energy will be split between the two nearest FFT result bins (two different frequencies).  The phase of your frequency peak might be slightly ahead or behind the phase of the slightly lower FFT frequency bin. If so, it will have the opposite phase relationship (slightly behind or ahead) of the the FFT frequency bin slightly on the higher or opposite side in frequency difference.

So it's not two different phases.  It's two opposing phase relationships when your dominant frequency spectra is correlated against the two nearest frequency FFT basis vectors, that are both slightly different (and one integer period of the FFT length apart), but on opposite sides in frequency, and thus opposite in phase difference relationship to your dominant frequency peak.

Unless, of course, your dominant peak represents a signal that is exactly integer periodic in the FFT aperture, in which case only one FFT result bin (per harmonic) will be involved, and you won't see two phases (for any non-zero FFT result bins that are above the noise floor)


Additive White Gaussian Noise (AWGN) and Undecimated DWT

One property of Orthogonal Transformations is that White Noise stays White (Uncorrelated) under Orthogonal Transformations (One could say it's a property of White Noise).

Usually this property assist with deriving properties of operators (Things are simpler).

Yet this property doesn't not hold for Transformations in general which might "Color" the noise in the new base.


Examples of linear/logistic regression applications

In image processing all I can think of is using linear regression to fit a plane to the intensity surface of an image, and using logistic regression for some  kind of supervised thresholding...

You might be better off using an example from computer vision instead. For example texture classification or object recognition.


How can I load a wav file, split it in frames, apply different filters do different frames and then recombine the filtered frames into a single file?

add technique using ann "hann" window and it works quite well but between the different frames I'm not able to avoid the discontinuities

That's probably because the filtering has tails that extend beyond the frame.  The window is supposed to taper the signal to 0 at the edges, but after filtering it is no longer 0 at the edges, which introduces discontinuities.

Are you zero-padding before filtering to keep these response tails?

I believe you need to apply a sqrt(window) before filtering ("analysis window"), and then another sqrt(window) after filtering ("synthesis window"), so the combination is a normal window, and the tails are tapered to 0, and then overlap the windows by the appropriate amount so that they sum back to unity.  Maybe also need to zero-pad before filtering so that the filtering tails don't wrap around, if you're doing FFT filtering, for instance.  More information here: http://www.dsprelated.com/dspbooks/sasp/Weighted_Overlap_Add.html

(If you're doing positive-phase time-domain filtering, it might be sufficient to use a regular window, zero-pad before filtering, and then instead of windowing again, just keep the entire zero-pad region with the tail in it when you sum back together, but I haven't tried it so I'm not sure.)


That definition is rubbish.

An eigenfunction of a system is a function that passes that system, only scaled by a scalar value.

Now, suppose your system is a simple complex band pass filter, which has a pass band at $f_0$, but stop bands everywhere else, especially at $-f_0$.

Now your cosine $\cos(2\pi f_0 t)$ has spectral components at both $f_0$ and $-f_0$. The output of that filter contains only a multiple if the positive part. Hence, the output is not a scalar multiple of the input cosine. Done.


Is my solution to the DTFT of a delta function correct?

You are absolutely correct.  For more cool, calculate the DTFT of $\delta \left(n - k \right)$.


What are the steps for applying a 50% Overlap to an FFT?

For a 50% overlap, offset successive 4096 sample windows by 2048.  e.g. the first half of the data in each FFT will be duplicated from the last half of the previous FFT window (before applying the Von Hann windowing function to either).


Inverse of the Gibbs Phenomenon to generate a sine wave

Square waves are composed of sine waves at the odd harmonics- i.e. 1st, 3rd, 5th, etc.  So yes, you can produce a sine wave from a square wave if you filter it with a low-pass filter whose bandwidth is at least as big as the first harmonic frequency, and whose cutoff frequency is lower than the third harmonic.


Sinusoids are power signals, and in the case of power signals, cross-correlation between two signals $x(t)$ and $y(t)$ is defined by

$$R_{xy}(t)=\lim_{T\rightarrow\infty}\frac{1}{2T}\int_{-T}^{T}x(\tau)y(t+\tau)d\tau\tag{1}$$

If $x(t)=\cos(\omega_1 t)$ and $y(t)=\cos(\omega_2 t)$ then for $\omega_1\neq\omega_2$ the cross-correlation at $t=0$ is

$$\begin{align}R_{xy}(0)&=\lim_{T\rightarrow\infty}\frac{1}{2T}\int_{-T}^{T}\cos(\omega_1 t)\cos(\omega_2 t)dt\\
&=\lim_{T\rightarrow\infty}\frac{1}{2T}\left[\frac{\sin((\omega_1-\omega_2)T)}{\omega_1-\omega_2}+\frac{\sin((\omega_1+\omega_2)T)}{\omega_1+\omega_2}\right]\\&=0
\end{align}$$

Note that the definition (1) avoids all convergence problems that you would encounter otherwise.


Effects of averaged periodogram (Welch's Method) on resolution?

Yes, if you reduce the FFT window from 4096 to 1024 samples, then you will have 4 times less frequency resolution. Also, if you have significant overlap between the windows, your resolution will not change. It only depends on the length of the FFT window. Furthermore, with Welch's method, a window for each segment is usually used, which again reduces the spectral resolution due to convolution of the signal spectrum with the window spectrum.

So, yes Welch's method reduces noise, but also spectral resolution.

If you want to have both, why dont you increase the recorded signal length to lets say 4*4096? This would amount to 360ms, which, depending on your requirements, might still be fine for guitar pitch detection. THere's always the tradeoff between spectral resolution and required amount of time.


Strong (0,0) pixel among an otherwise normal cross correlation peak

Figures that the same day I post a bounty I answer my own question.

The answer to this question is dead/speckle like pixels. For a fixed number of dead pixels of a given brightness, the smaller the relative shift between the two images, the fewer dead pixels that are needed to cause the (0,0) pixel to spike. It can take fewer brighter dead pixels to have the same effect as more dimmer ones. Lastly, the threshold between a strong (0,0) pixel and normal one is extremely tight in the number of dead pixels. Literally a change in just 10 dead pixels can massively increase the probability that the (0,0) pixel becomes a problem.

This problem can be solved by a combination of dead pixel finding/interpolation and applying a light Gaussian blur to the images to minimize speckle.


Algorithms for computing FFT in parallel

I think your main problem isn't how to parallel the algorithm (which actually can be done) but it's numeric precision.
FFTs of a size that large are numerically quite tricky. The FFT coefficients are of the form
$e^{-j\cdot 2\cdot \pi \cdot \frac{k }{N}}$ and if N is very large the coefficient calculation gets noisy. Lets say you have $N = 2^{40}$ and you use 64 bit double precision arithmetic. The first 1000 coefficients have a real part that's exactly unity (although it shouldn't be that way) so you will need higher precision math, which is very inefficient and cumbersome to use.

You will also amass lots of rounding and truncation errors as the sheer number of operations that go into a single output number is very large too. Due to the "every output depends on every input" nature of the FFT, error propagation is rampant.

I'm not aware of an easy way to work around that. Your request is an unusual one. Most applications that do spectral analysis of large data sets do a running analysis where you don't have that problem. Perhaps if you can describe your application and it's constraints a but more, we can point you to a more suitable solution.


What is the intuition behind convolution?

No doubt there are many links and expositions on wikipedia and elsewhere about convolution, its relation to linear systems, and its link to frequency analysis, however I do not think that is what you are after.

Forget convolution for a moment. Just think of nothing but averages:

If I give you a series of say, $10$ numbers, and ask you to take the average, you can do that very easily.

Let us say that I give you another set of $10$ numbers again, and also ask you to take the average. You can again, do that very easily.

I keep giving you a set of $10$ numbers, and you keep giving me averages. What a life!

Now imagine I tell you a secret: Those sets of $10$ numbers I have been giving you, are actually from a taken from a long series of numbers. In fact I have a list of $100$ numbers. The first set of $10$ numbers I gave you were the first $10$ numbers of that list. (First sample to 10th sample). You averaged them.

The second list of $10$ numbers were actually the 2nd sample, to the 11th sample. You averaged them.

The third list of $10$ numbers I gave you were actually the 3rd sample, to the 12th sample. You averaged them.

After all this, I ask you to plot all the averages you gave me, one after the other.

Congratulations! You just showed me a convolution result!

(Disclaimer: There are some border conditions that I have skipped over, but this is basically all a convolution is).

So in summary, think of a convolution as just a sliding average.

Now, remember how I said "average" those $10$ numbers I gave you? Well you can also do a weighted average. (That is, instead of multiplying all the $10$ numbers by $1$ before averaging them, we multiply them by different numbers, before averaging them. This is a weighted average). This will still be a convolution though. We call those numbers "weights", and sometimes they are referred to as the "filter kernel".

Some advanced systems actually change their weights as they run, and they are called adaptive filters.

Eitherway, this should provide you with a solid intuition: Just remember that convolution can be thought of - is - a moving average of your data. Hope that helped!


Compress a signal by storing signal diff instead of actual samples - is there such a thing?

You can also think of delta encoding as linear predictive coding (LPC) where only the prediction residual ($x[n]-\hat{x}[n]$ in @robertbristow-johnson's notation) is stored and the predictor of the current sample is the previous sample. This is a fixed linear predictor (not with arbitrary coefficients optimized to data) that can exactly predict constant signals. Run the same linear predictive coding again on the residuals, and you have exactly predicted linear signals. Next round, quadratic signals. Or run a higher-order fixed predictor once to do the same.

Such fixed predictors are listed in Tony Robinson's SHORTEN technical report, yours in Eq. 4, and are also included in the FLAC lossless audio codec although not often used. Calculating the best prediction coefficients for each data block and storing them in a header of the compressed block results in better compression than the use of fixed predictors.

For $m$-bit input the residual is an $m+1$ -bit number, because it is the difference of an $m$-bit input and an $m$-bit prediction. However, removing the most significant bit (MSB) of the residual has no consequence in $m$-bit modular arithmetic, so the residuals can be stored as $m$-bit numbers.

The linear predictor is supposed to do the whitening, making the residuals independent. In lossless compression, what is left to do is to entropy code the residuals, instead of using run-length or other symbol-based encoding that doesn't work so well on noisy signals. Typically, entropy coding is done by a prefix code (also known as prefix-free code) that assigns longer code words to large residuals, approximately minimizing the mean encoding length for an assumed distribution of the residual values. A Rice code (also known as Golomb–Rice code or GR code) variant compatible with signed numbers can be used, as is done in FLAC (Table 1), or signed exp-Golomb code as is done in the h.264 video compression standard. Rice code has a distribution parameter that needs to be optimized for the data block and saved in the block header.

Table 1. Binary codewords of 4-bit signed integers encoded in Rice code with different Rice code parameter $p$ values, using FLAC__bitwriter_write_rice_signed (source code). This variant of Rice code is a bit wasteful in the sense that not all binary strings are recognized as a codeword.
$\begin{array}{rl}
\begin{array}{r}\\-8\\-7\\-6\\-5\\-4\\-3\\-2\\-1\\0\\1\\2\\3\\4\\5\\6\\7\end{array}&\begin{array}{lllll}
p=0&p=1&p=2&p=3\\
000000000000001&000000010&000110&01110\\
0000000000001&00000010&000100&01100\\
00000000001&0000010&00110&01010\\
000000001&000010&00100&01000\\
0000001&00010&0110&1110\\
00001&0010&0100&1100\\
001&010&110&1010\\
1&10&100&1000\\
01&11&101&1001\\
0001&011&111&1011\\
000001&0011&0101&1101\\
00000001&00011&0111&1111\\
0000000001&000011&00101&01001\\
000000000001&0000011&00111&01011\\
00000000000001&00000011&000101&01101\\
0000000000000001&000000011&000111&01111\end{array}\end{array}$

As a further enhancement, encoding not just one but multiple residuals into a single codeword can more accurately accommodate the true distribution of residuals and may give a better compression ratio, see asymmetric numeral systems.


response of LTI system to a Random Input Signal

There are some basics that you really need to understand. The expectation operator of a deterministic (non-random) function just equals the function. Since the impulse response $h(t)$ is considered to be deterministic, you have

$$E[h(\tau)X(t-\tau)]=h(\tau)E[X(t-\tau)]$$

Furthermore, since $X(t)$ is known to be stationary, its mean doesn't depend on time, so you have

$$E[X(t-\tau)]=E[X(t)]=m_x$$


Define this

  Tn=0:Ts:1;

to be

  Tn=0:Ts:1 -Ts;

By defining it the way you have in the program you have not included an integer number of cycles in the samples. This will result in spectral leakage because your sinuoid now does not fit into one frqeuency bin and thus spill out into other bins and this reduces in magnitude.


Significance of $ \lambda $ in Basis Pursuit

There are 2 forms of the Basis Pursuit problem:

$$\begin{align*}
\text{The $ \lambda $ Form:} & \quad && \arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} + \lambda {\left\| x \right\|}_{1} \\
\text{The $ \epsilon$ Form:} && \quad & \arg \min_{x} {\left\| x \right\|}_{1} \\
& \text{subject to} &&  \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} \leq \epsilon
\end{align*}$$

Sometimes it is more intuitive to have a look on the $ \epsilon $ form as it tell us exactly the trade off.
The lower $ \epsilon $ it means we require teh solution to hold the Linear Equation and we pay for that with the $ {L}_{1} $ Norm of the solution.

The nice things is that for any value $ \epsilon $ we can find $ \lambda $ for the other form in such way the generate the same solution:

$$ \forall \epsilon, \; \exists \lambda : x \left( \epsilon \right) = x \left
( \lambda \right) $$

Where $ x \left( \epsilon \right) $ and $ x \left ( \lambda \right) $ are the solution of the $ \epsilon $ form and the $ \lambda $ form respectively.

You may see an example in my answer at Cross Validated Question 291962.


Is it possible to combine two or more Biquads?

Possible, yes.  But higher-order IIR filters tend to be less stable and more susceptible to various numerical problems (lost precision, underflow, overflow, etc.), compared to a cascade or bank of suitable partitioned biquad filters.


Why do we need to minimize the TX power to require for meeting QoS?

I interpret that sentence the following:

  Under energy harvesting contsraints and Quality of Service constraints, the second category of works strive to minimize TX power.

So, TX power isn't minimized in order to meet QoS constraints, but as far as possible while still meeting QoS constraints.


This sounds like a confusion in terminology. See http://en.wikipedia.org/wiki/Digitizing

Digitization basically involves two steps:

Discretization: Sampling the signal at discrete times
Quantization: Turning the samples from (in theory) infinite resolution to finite resolution

Look at all possible combinations you can have 4 different types of signal:

continuous time, continous amplitude
discrete time, continous amplitude
continuous time, discrete amplitude
discrete time, discrete amplitude

Typically we refer to #1 as "analog" and #4 as "digital". #2 and #3 are of little practical relevance, other than getting the mathematical models right. Quantization and Discretization are really quite different mathematically so it doesn't make a lot of sense to analyze them in combination. In other words: the theory of time discretization applies equally to continuous amplitude and discrete amplitude signals.


Determine fundamental period from zero crossing signal

In general, you can not determine period of a signal from just the zero crossings.  If the signal has harmonic content, then a pitch detection/estimation algorithm, such as autocorrelation, might be one solution.  If the signal is more spectrally pure (e.g. close to sinusoidal fundamental, little overtone or harmonic energy), then a suitable DSP band-pass filter applied before locating the zero-crossings may help get rid of "noisy" samples.


Using Least Squares for curve fitting in Matlab

1) Form a matrix, A of order $m$ x $3$, such that:

$$A=\left[f(x) \ \ g(x) \ \ h(x)\right]$$

where $f(x), g(x)$ and $h(x)$ are three column vectors such that each column is a sample of component guassian at $x$.

2) Form a column vector, $b$ of your original Gaussian pulse of order m.

3) Solve the system $x=((A'*A)^{-1}*A')*b$, where $x$ is your desired set of coefficients of order $3$ x $1$, that would minimize the error in the least square sense.


IMHO the best book that relates circuits to the z-transform is Papoulis, "Circuits and Systems: A Modern Approach". You can find it used for a handful of dollars.


The question as stated doesn't make sense. In particular, "sixteen frequencies per bit time" doesn't mean anything. So, I'll try to guess at what the question actually is, and try to point you in the right direction.

Let's break the question in two parts. First, we're told that the carrier can have eight different amplitudes. This means we can encode three bits in the carrier amplitude. If we label amplitudes $L_0,\ldots,L_7$, we can encode the bit sequence $000$ as amplitude $L_0$, $001$ as amplitude $L_1$, etc. This means that this carrier can transmit bits 3 times faster than a carrier with only two allowed amplitudes.

The second part of the question seems to indicate that you actually have sixteen carriers, all operating simultaneously. (and that they don't interfere with each other). Each carrier can have eight different amplitudes.

If you take $R$ as the baseline rate obtained by a single carrier with two amplitudes, you should now be able to calculate the rate of 16 carriers with 8 amplitudes each, as a multiple of $R$.

Edit: I think there is a second possible interpretation of the question: this could be some sort of combination ASK-FSK, where there are 16 possible frequencies (FSK) and each could have 8 possible amplitudes (ASK). Each frequency would carry 4 bits and each amplitude a further 3 bits, for a total of seven bits per symbol time.


In voice recognition system, how to match digital input to a given word without seconds delay?

There is a massive range of possible sounds that a human might recognize (or mis-recognize) as the same word in their language.  The problem is searching a large dimensional space representing all those possibilities (continuous variations of spectral and timing differences) for the closest match, perhaps also given some context around the search target for the most probable of multiple near matches.  This is a computational complexity problem that can sometimes be reduced in time by more computational power (a 1000 processor server or supercomputer might be able to do it in less than "seconds" of delay).

Or you can search over less degrees-of-freedom with less robustness or accuracy, and perhaps miss most possible word matches.


Continuous and discrete-time raised cosine filter properties in frequency domain

Root Raised Cosine has a bandwidth $B=\frac{1+\alpha}{2T_c}$ which is $(1+\alpha)$ times more than the minimum required bandwidth set by Nyquist no-ISI criterion. To satisfy the Nyquist sampling theorem now, one must sample at a rate of $2B$. In time domain, we can say that we have to have more samples/symbol than 1. In practice, and to facilitate the synchronization blocks down the chain, this sampling rate is selected such that we have an integer number of samples/symbol. The smallest integer greater than 1 is 2. In your eq (6), you are sampling at a rate of 1 sample/symbol which induces aliasing. Hence, the results of eq (5) and eq (6) will not be equivalent. For 2 or more samples/symbol for incoming signal and assuming no analog distortions, eq (5) and eq (6) will become equivalent after downsampling to symbol rate.

The best way to understand is to imagine a figure of RRC in frequency domain with aliases at integer multiples of symbol rate $1/T_c$.


Is autospectrum the same as power spectrum density?

This may be just a matter of semantics. Indeed it is true that "The Fourier Transform of the auto correlation x[n] is identical to the magnitude squared of the Fourier Transform of x[n]". This is just a mathematical identity.

You could call that the power spectrum density (PSD) but in most practical applications any actual PSD computation will involve some sort of framing and windowing. The choice of framing and windowing will impact the result (and it's a complicated trade off), so there isn't really one clear and unambiguous definition of PSD.

You can still use the mathematical identity but it needs to be properly adjust with respect to framing, windowing and circular vs. linear convolution.


Estimate/Compute parameters of image dataset

Height is easy. You know that the camera is at $(0,0,0)$. Just fit a plane to the ground pixels and measure the orthogonal distance from the origin to this plane.

Frame rate is more of an issue. It is correlated with the speed of the user moving the camera. I suppose an approach would be to actually take the Kinect camera, move it like a normal person does and record the frame rate. You could then compute the optical flow in your dataset and record the magnitudes of flow vectors. You could then train a classifier (e.g. random forest) to learn to regress the frame rate given the flow vectors. Assuming that the camera motion is very similar, you could then put your dataset into this classifier and regress the frame rates for the whole trajectory. The average over all sets of frames should give you, kind of, the frame rate.


Starting with a Bandpass signal - this is a signal that has a bandwidth $BW$ about a center frequency $f_c$, so the signal content goes from $f_c-1/2BW$ to $f_c+1/2BW$. For a real signal (i.e. not complex) the their is also the corresponding negative frequency component i.e. the signal also exists at $-f_c$.

Basebanding often refers to shifting this signal down to DC. This is done by mixing the signal (modulating) it with $\exp(-j2\pi f_c t)$ and low pass filtering to remove the component that was originally at $-f_c$ and gets moved to $-2f_c$ due to the shifting. The component that was at $f_c$ and gets moved to DC is kept by the low pass filter. This signal is now complex due to the loss of conjugate symmetry. This resulting signal is referred to as the Baseband signal.

Note that the mixing signal $\exp(-j2\pi f_c t) = \cos(j2\pi f_c t)-j\sin(j2\pi f_c t)$, which is often how you see it implemented in communication systems.

You can achieve the same result by using the Hilbert transform which essentially removes the negative frequency component and then shifts the signal down to baseband.

References:
Haykin - Communication Systems (but almost any communications textbook), also Crochiere - Multirate Digital Signal Processing (somewhere in the first couple of chapters they explain this).


How can I plot Wiener filter's behaviour for different taps (memory filter)?

Only action needed is to change the value of M and study the plot .


Accuracy of Goertzel algorithm for measure harmonic

The most common Goertzel algorithms are identical to calculating the magnitude of 1 bin of an FFT/DFT using a rectangular window.  As such, there will be rectangular window artifacts (so called "leakage") for any spectral content, such as harmonics, that are not exactly periodic in the DFT length, or Goertzel length.

Thus the accuracy of using just 1 bin may not be more accurate if you increase the length from one where all the spectral content is periodic to one where it isn't.  With a complete FFT, one can somewhat mitigate this effect by using high quality interpolation in the complex domain before taking the magnitude.


Cannot recreate a simple Nyquist series in Matlab

HINT:

Try Y1 = cos(2*pi*fNy*X);

(And think about the $x$-axis values at which you sampled the sine)


Generating a Sampled Square Wave in Matlab

The word "points" means "samples." Product fs*T is fs samples/second times T seconds. That's a signal duration of 40*5 = 200 samples (or points). This means your Matlab code for vector t should have been:

t = [0:1/fs:T-1/fs];

You wrote tau = 1/100, so (1–tau) = 0.99. Ratio fsT(1-tau)/2 = 40*5*0.99/2 = 99 samples (points). You cannot have a square wave sequence of length 200 samples (points) if the first and last 99 samples (points) are zero-valued. Your fs, T, and tau variables are incompatible with each other!

Note: Your x = square(2*pi*t); Matlab command is equivalent to:

x = square(2*pi*1*t); % Note the "1"

which produces square wave that repeats every one second. If you want just one cycle of a square wave over 200 samples use the commend:

x = square(2*pi*t/5);


The random demodulator (in this definition) compresses the signal by taking samples that are each linear combinations of W/R Nyquist rate samples. So it does not take W/R samples at a time; it takes one sample at a time, each of which contain information corresponding to several (W/R) Nyquist rate samples. It is then up to a reconstruction algorithm to subsequently "disentangle" these into the original signal (approximately) at Nyquist rate.

The random demodulator works by "smearing" the content of the input signal all over the spectrum (by multiplication with the spreading sequence). Provided that the input signal is sufficiently sparse, we can now reconstruct it from just a small portion of the resulting spectrum (selected by the low-pass filter).


It depends on your decoder; for example, in message-passing decoders, most commonly, a prefactor doesn't affect the results at all, so you can scale to your heart's delight.
It's also common to just clip – whether a bit has a LLR of -100 or -90 doesn't really matter, does it.

However, the actual confidence levels don't make sense then anymore.

Well, you'd have to transform them, too. Notice that scaling LLRs makes sense in a Gaussian noise setting, as LLR involves taking the logarithm of a $ e^{-a|x-x_i|^2}$ function, so you're left with something rather sensible (scaling the LLR just scales $a$, which is common to all samples and constellation points). I don't know how your noise looks or how you arrive at confidence levels, but try to apply the logarithm to a ratio of confidences and see whether scaling the LLR preserves sensibility.


No you are not missing anything. A moving average with a period of $T = 1/60$ seconds will indeed have a notch at 60 Hz, since the frequency response is a Sinc function with the first null at 1/T.  (So it is both a low pass filter with an magnitude envelope that goes down as 1/f combined with a notches at integers of 1/T not including 0).

This is because a moving average has an impulse response that is a rectangular function of duration T, and the Fourier Transform of such a function is a Sinc (and the Fourier Transform of the Impulse Response is the Frequency Response). Digital Filters are similarly formed when the impulse response is a sampled rectangular function and thus their Fourier Transform will approach a Sinc function (what I like to refer as an "aliased" Sinc function similar to other aliasing in discrete time systems). If we position the sampling rate at integer multiples of 60 Hz, starting at 120 Hz and higher, the null will be at 60 Hz and it's higher harmonics if the digital impulse response (formally called the unit sample response) is uniform samples for duration T, where T is 1/60 Hz (which is a moving average specifically when the weighting of each sample is 1/N for N samples, otherwise it is a scaled moving average). Therefore the minimum sampling rate required to do this digitally would be 120 Hz, otherwise 180 Hz, 240 Hz, and higher multiples of 60 Hz can also be used).  Note that I said this would approach a Sinc function; notice the plots below of the actual frequency responses for each choice of sampling rate when we position the null at 60 Hz and do the moving average, given by the actual digital transfer function according to:

$$H(z) = \frac{1}{N}\sum_{n=0}^{N-1}z^{-n}$$

With the frequency responses below determined by using the unit circle for z ($z= e^{j\omega}) $ where $\omega$ is the normalized angular frequency from 0 to $2\pi$ corresponding to 0 to the sampling rate. The actual Matlab/ octave command that can be used to readily determine this is as follows (and Python SciPy has a similar command):

freqz([ones(N,1), 1/N]);

Where N for each of the cases below is 2, 3, 4, and 10.

Also what you could alternatively consider which is even simpler is a "comb filter" structure which is simply the addition of the input with a delayed version of itself where the delay is equal to T/2.

This will have the frequency response according to:

$$H(\omega) = 1 + e^{-j\omega T/2}$$

Since the Fourier Transform of a delay is given as

$$x(t) = \delta(t-T) <==> x(\omega) = e^{-j\omega T}$$

With the magnitude plotted below:

This is equivalent to the digital filter with a delay of one sample when the sampling rate is 120 Hz, and a transfer function given as:

$$H(z) = 1+z^{-1}$$

with as stated previously the frequency response when z is the unit circle as

$$H(e^{j\omega}) = 1 + e^{-j\omega}$$

with $\omega$ going from 0 to $2\pi$ in this case representing the normalized angular freq from 0 to the sampling rate.

Like the moving average case above, this could also be implemented at any sampling frequency that is an integer multiple of 60 Hz, starting at the minimum 120 Hz rate.

If you are concerned about the "frequency droop" in your passband area for either of the above two solutions, you could consider using this very simple three tap FIR compensator that I described at this link to help flatten the passband:

how to make CIC compensation filter

In all the cases above, the sampling rate must be an integer multiple of the notch frequency desired. If that is too restrictive, here is a digital notch filter that I would typically implement when a very tight notch is required with a flat passband, to compare to the approaches above (and does not have the same limitations on sampling rate):

Transfer function of second order notch filter

Side note: The relatively high sidelobes, or stated another way, the slow frequency roll-off of 1/f in magnitude, is why a moving average is often a poor choice as an FIR filter for low pass filtering (or estimating the DC value which IS the average) UNLESS you are filtering what is known to be white noise. For the case of white noise it is one of the best choices. Why? Because the moving average does have the narrowest main lobe for the given number of taps compared to any other FIR filter, and if the noise is known to be white, then the 1/f magnitude roll-off will be sufficient to suppress the higher frequency components.


That's a really hard problem. If you have access to the original snare stem and the amount of post processing (effects, panning, mastering) is moderate you may be able to pull some of it out with an adaptive filter.

There are a few commercial plug ins that do this type of thing, but they all have their strengths and weaknesses so it depends a bit on your specific needs


Decoding of binary signal from a temperature sensor

It is impossible to extract the correct data from this string with only one sample (as is mentioned by Laurent Duval in the comments). Here is something that caught my eye:

Markers

You have these four blocks at the beginning and the end of the signal that are some kind of start/end markers, I suppose. Are you sure that you caught the whole markers?

1000111110010010110001100011111111000 could be split up in the following way:

1000
1111

100100101100011000111

1111
1000

so the 1000 and 1111 might actually be part of the markers.

The temperature itself

You say the signal corresponds to 21.0°, i.e. a number with an accuracy of 0.1°. It could be that this is encoded as the integer 210. 210 in binary is 11010010 and you can find this string (reversed) in there:

100100101100011000111
  01001011

The prefixing 10 is 2 in decimal and could be a 2 bit indicator of the number of 4 bit words that encode the temperature, however this would be strange since the temperature is reversed and the prefix is not...

Please try the following:

Measure the same temperature multiple times in independent experiments. Is the resulting data string the same in each experiment? Are parts of the string the same in each experiment?
Measure multiple temperatures that only differ by 0.1°. Is there a corresponding pattern in the data string?
Do you have a datasheet of your sensor?


What is the advantage of Wiener filter for noise reduction of a 1D signal?

The Wiener filter considers statistical behaviours of the noise and the signal, and thus, (theoretically) achieves optimum separation of them for a class of signals and systems, which is not the case for more classical approaches.
Wiener filter frequency response is such that, at those frequencies where noise power is dominant (a.k.a. low SNR), the gain is reduced, and the output is suppressed; causing noise-reduction. Whereas for frequencies at which the signal power is dominant (a.k.a. high SNR), then the gain is closer to one, and output is closer to the input.
Thus, to determine the gain of the Wiener filter, at a particular frequency, you must know the Power Spectral Densities (PSDs) of the noise and the desired signal.
Typically, you won't know these in advance very well; but they can be estimated from available data. The result is a departure from the optimum (ideal) performance.


What could cause fast Fourier transform to give complex conjugate of the intended result?

Taking the inverse FFT (IFFT) instead of FFT would cause a phase inversion.


You have to carefully read the documentation of the functions that use/generate normalized frequency in some cases they use $fs$ and in others they use $f_s/2$.
The Mathworks docuementation for firpmord:
"Description
[n,fo,ao,w] = firpmord(f,a,dev) returns the approximate order n, normalized frequency band edges fo, frequency band amplitudes ao, and weights w that meet input specifications f, a, and dev.
example
[] = firpmord(,fs) specifies a sampling frequency fs. fs defaults to 2 Hz, implying a Nyquist frequency of 1 Hz. You can specify band edges scaled to a particular application's sample rate. You can use this with any of the previous input syntaxes."
In this case the maximum frequency that can be given in the f vector is 1.0 which corresponds to $f_s$. This is likely because the Matlab code is designing real (not complex) filters and you can only specify the response from 0 -> $f_s/2$ since the other half will be determined by conjugate symmetry.
It can be confusing - looking at Matlab's Remez and Remezord functions, they tend to say they normalize by the Nyquist frequency (they mean $f_s/2$) but they don't explicitly tell you that. Looking at the graphs of the filter responses in their examples in the documentation does make it clear - they show only half the frequency range and the maximum frequency is 1.0 (meaning it has been normalized by $f_s/2$.
In other code, e.g. ones that deal with complex filters, you may see the normalization use $f_s$, because you no longer have the conjugate symmetry you can specify the response over the whole frequency range (0-> $f_s$).
So it really depends. There are really two conventions for normalized frequency and you have to read the documentation or ask what whether they've normalized by $f_s$ or $f_s/2$.


Proof of properties of Fourier series in CT

Looking at the change of terms in the equality:
$$\frac{a}{T}\int_{T/a}x(at)e^{jk(\omega_0a)t}dt = \frac{a}{T}\int_{T}x(L)e^{jk\omega_0L}\frac{dL}{a}$$
$$at = L\implies \begin{cases}\text{if }t=T/a\implies L=T\\adt=dL\implies dt=\frac{dL}{a}\end{cases}$$
Looking at the equality again:
$$ \sum^{+\infty}_{k\ =\ -\infty} \sum^{+\infty}_{l\ =\ -\infty} a_k b_le^{j(k+l)\omega_0t}=\sum^{+\infty}_{m\ =\ -\infty} \left[\sum^{+\infty}_{l\ =\ -\infty} a_{m-l} b_l\right]e^{jm\omega_0t} $$
$$m=k+l\implies \begin{cases}k=m-l\quad({\scriptstyle{\rm indeed}})\\\text{and } m-l=\pm \infty \implies  m=\pm\infty+l=\pm\infty\end{cases}$$


What is the difference between wide sense and strict sense stationary processes?

A process is stationary if:

its mean is a constant value: $\mu_x(t)=\mu x$
its MSV(mean square value) is a constant value.
its variance is a constant value. $\sigma^2_x(t)=\sigma^2_x$
its autocorrelation depends on the time distance between 2 samples : $$R_x(t_1,t_2)=R_x(\tau)$$
its auto-covariance depends on the time distance between 2 samples at time $t_1$ and $t_2$ : $$K_x(t_1,t_2)=K_x(\tau)$$
So, with respect to information above, we cannot just say that a process is stationary if it does not depend on time, because for mean,MSV and variance a stationary process does not depend on time since its values are constant, but for autocorrelation and auto-covariance it does depend on time actually, but the function itself does not depend on time since it depends only on the time distance between 2 samples.

The process is SSS (strict-sense stationary) if its $N^{\rm th}$ order probability density function is stationary for any given $N$.

The process is WSS (wide/weak-sense-stationary) if its mean value is constant $\mu_x(t)=\mu_x$, and its autocorrelation function depends on time difference between 2 samples, $$R_x(t_1,t_2)=R_x(\tau).$$

Important:
Remember that all SSS processes are also WSS, but NOT all WSS processes are SSS.There is only one exception: "Jointly Gaussian Processes". Jointly Gaussian processes have any order of probability density functions(property of SSS). So, even if a jointly gaussian is WSS it is also SSS.And also remember that white gaussian noise process is a jointly gaussian process with 0 mean. So it is WSS and also SSS. Since its mean value is 0, the we can say that its autocorrelation function and autocovariance function are equal and they do depend on time difference between 2 samples  at time $t_1$ and $t_2$. $R(\tau)=K(\tau)$ where $\tau=t_1-t_2$.

In conclusion,

if you see only the mean value is constant (not depend on time), and the autocorrelation function depends on only time difference $\tau= t_1-t_1$. $R(\tau)$. --> WSS
if you see that the $N$-th order pdf of process is stationary for any $N$ --> SSS
if the process is jointly gaussian --> WSS and SSS
if the process is white gaussian noise process --> WSS and SSS with mean=0 and $$R(\tau)=K(\tau).$$


How to Compensate for Frequency Offset in single carrier transmission Using Coarse and Fine Compensation?

To compensate for a frequency offset first it needs to be measured, let us take the example of  PSK.

At the trasnmitter : Suppose the signal is upconverted at the transmitter to $F_c$, then the received signal is given by $$s(t)e^{j2\pi F_ct}$$ here $s(t)$ is the baseband PSK signal.

At the Receiver: Suppose the received signal has a frequency offset, either due to doppler shift or inaccuracy of the receiver clock, let this offset be $\Delta_f$, after downconversion at the receiver, the signal sent to the ADC is $$s(t)e^{j2\pi F_ct}s^{*}(t)e^{j2\pi( F_c - \Delta_f)t}$$

$$\vert s(t) \vert^2 e^{-j2\pi\Delta_ft}$$

Where for the sake of simplicity, let us normalize $$\vert s(t) \vert^2 = 1$$

After going through the ADC the digitized samples are given by $$ y[n] = e^{j2\pi\Delta_fnT_s}$$ where $T_s$ is the sampling interval

Let us now correlate adjacent samples of this digitized signal, we get

$$ y[n+1]y[n]^* = e^{j2\pi\Delta_f(n+1)T_s}e^{-j2\pi\Delta_fnT_s}$$

$$ y[n+1]y[n]^* = e^{j2\pi\Delta_fT_s}$$

Everything other than $\Delta_f$ is known in the above equation, which can be found out using the following relation, let us denote the phase resulting from the correlation of samples as following:

$$ y[n+1]y[n]^* = u_k$$ then
 $$ \Delta_f = \frac{\Phi_k}{2\pi T_s} \tag{1}$$

Where $$\Phi_k = arctan\left( \frac{Img(u_k)}{Real(u_k)}  \right)$$

For a good result which is less impacted by noise, take average over a larger number of such correlations between incoming samples.

For the case when $\vert s(t) \vert^2$ is not unity or unit power, the value $E(s[n+1]s[n])$ is decided by the constellation. Usually a pilot is sent with a known constellation point to the receiver and trasnmitter to determine the frequency offset and channel.

For ex: in 4 QAM, you could send the constellation point/symbol $\sqrt{\frac{1}{2}} + j\sqrt{\frac{1}{2}}$

If however the SINR is low, the envelope of the signal $\vert s(t) \vert^2$, might not be predictable at the receiver in such cases it is better to work with the demodulated symbols and not the samples. The upside is better offset estimation, the downside is a longer convergence loop, because of need to demodulate and then  recover offset.

In such cases, the receiver will first try to demodulate the symbol with the frequency offset (hopefully it's not too big and can be demodulated correctly). The rotation in complex plane is directly linked to the offset as in equations above.

Now $T_s$ would be replaced by the symbol time in the above equation (1), because we would be correlating between adjacent symbols and not samples, rest of the analysis will remain exactly the same.

Offset correction at the receiver:

Once the frequency offset is estimated it can then be corrected at the receiver by applying the adjustment to the PLL at the receiver.

However, the decision to correct the offset also depends on the SINR and quality of measurements. If the SINR is low (that means the measured frequency offset is not fully trust worthy, then correct only a part of the measured frequency offset, ex : 75% in one  estimation cycle, and successful estimations will then converge to the actual offset over a period of time). So the amount of total correction once the estimation of the offset is done should be reflected in the confidence in the estimate decided by the SINR and other physical layer parameters.


FFT on continuous fifo data - Is y axis magnitude sum of all samples?

In general the  output of a Fourier Transform is "spectral density". For a continuous Fourier Transform that's implicit: if the signal in Volts, the unit of the FT are Volts/Hz. For the discrete case, the units are the same but the frequency interval is implicit: it's the sample rate divided by the FFT length.
If you want to extract "physical" meaning from the FFT, you need to watch your scaling. In your specific case it would probably be best to scale both the forward and backward FFT by $1/\sqrt{N}$ since that maintains Parseval's theorem, i.e.
$$\sum |x[n]|^2 = \sum |X[k]|^2$$
which means that total power in the time and frequency domain are the same.


Matched filter and cross correlation in digital processing

I'll illustrate with an example: detecting a rectangular pulse in noise, in Matlab.

Let's start with defining the signal that contains the pulse.

fs = 1/1000;
t = 0:fs:1;
s = [zeros(1,400), ones(1,100), zeros(1,501)];

First let's do detection in noise with small power. We calculate the received (noisy) signal:

np = 0.1;  % noise power (variance)
n = sqrt(np)*randn(1,1001); % noise signal
r = s + n; % received signal
plot(t,r); title('Received signal, small noise power');

Now let's run the correlation:

c = xcorr(r,s);
plot(c); title('Correlation, small noise power');

You can see how the correlation produced a large output, indicating that the pulse was there. What is amazing is that the correlator can find the pulse even when it is impossible for one to see it. Let us increase the noise power and repeat the process.

np = 4;
n = sqrt(np)*randn(1,1001);
r1 = s + n;  % received signal with large noise

As you can see, the correlator has an incredible ability to detect signals that are embedded in strong noise.

One thing you can do for yourself is run the correlation when the pulse is not present. The signal will be hard to distinguish from the noisy pulse shown above, but you'll see that the correlation does not produce a sharp maximum.

So, the implementation is quite easy; what is not so easy to do is:

Prove that this technique produces the best possible signal to noise ratio at the correlation's peak.
Determine the optimum threshold given some desired constraints on the probability of wrong detection (false positive or false negative).


Typically you can improve performance by analyzing the data that you want and analyzing the data that you don't want. I would assume that a plane can't turn infinitely fast so the actual "want" data is properly fairly band limited. Depending on how fast you are sampling the accelerometer you may pick up a lot of "out of bandwidth" noise. This could easily be improved with a low pass filter or some moving average contraption.  The more you know about the properties of the actual signal and the properties of the noise, the better you can separate the two.


Transfer function estimation from frequency response

Here you find an extensive discussion about transfer function estimation and even source code.
Your problem can be expressed as $F(2j \pi f_i) \approx g_i$, where $f_i$, $g_i$ are your measurements, and you also you  can expresss $F(s) = N(s) / D(s)$ as a parametric function, then the parameters may be adjusted with curve fit.
Python provides curve_fit, that can be directly applied to this problem as in this answer, with the difference that you want to apply for x in the imaginary axis. This can be improved by also passing the gradient in terms of the parameters as well.
Maybe you prefer to express your function in terms of poles and zeros, this may be specially convenient if you want to ensure stability (you can add to your cost function the Lagrangian Multiplier for $real(p) < 0$)
def pztf(omega, h, p, z):
  num = np.prod(1j*omega[None,:] - p[:, None], axis=0);
  den = np.prod(1j*omega[None,:] - z[: None], axis=0);
  return h * num ./ den

Maybe you want to fit the logarithm of this, this may produce a better fit in the low gain frequencies.
There are also more analytic approaches. The so called Levy method is interesting,
$F(s) = \frac{N(s)}{D(s)}$ is equivalent to $F(s) D(s) - N(s) = 0$ this can be expressed as least square fitting.
This is biased in the sense that the errors close to the poles receives less weight, this can be mitigated by solving using $D_0(s) = 1$, and iteratively refining it by iteratively fitting
$$\frac{1}{D_{i-1}(s)}\left( F(s) D_i(s) - N_i(s) \right) = 0$$
when it converges after some iterations we have $D_{i-1}(s) \approx D_i(s)$ we have
$$\frac{1}{D_{i-1}(s)}\left( F(s) - N_i(s) \right) \approx F(s) - \frac{N_i(s)}{D_{i-1}(s)}$$
That corresponds to the original problem (without the $D(s)$ factor).
If you want better accuracy in different frequencies. You achieve this by multiplying the whole expression by a given $W(s)$, then te fitting iteration becomes
$$\frac{W(s)}{D_{i-1}(s)}\left( F(s) D_i(s) - N_i(s) \right) = 0$$


Is DCT (Discrete Cosine Transform) of Type-2 lossless or lossy?

Mostly yes, but it depends on the context. Let us elaborate.

DCT-II is one of the many forms of Discrete Cosine Transforms, and probably the  most widely used one, as it is (somehow) present in JPEG or MP3 formats. "Lossy" often refer to the compression standard which uses it, because the  main loss results from quantization (and generally not the transform by  itself). By itself, the transform is invertible, even more orthogonal or orthonormal, so in theory you have no loss (it is bijective).

As @Fat32 points out, its coefficients sometimes involve cosines, hence irrational numbers (except at specific values) that are not easily  represented in finite arithmetic (float or int), and can induce round-off errors.

But:

DCT-II of size $1$ is $1$, DCT-II of size 2 (unnormalized) is: \begin{smallmatrix} 1 & 1 \\1 & -1\end{smallmatrix} so DCT-II can be lossless (in special cases though),
For integer data (like images), it seems possible, with enough bit-depth computations (related to the initial integer range and the DCT size), to remain practically lossless, as you know that the  transform/inverse transform result should be integer,
Very  accurate integer, or dyadic approximations exist, meant for integer hardware: for instance the binDCT.

My final answer is thus: DCT-II is mostly lossless, if you take care of it.


pixels in the same image are correlated but how can I visualize it?

If you would like to show that in natural images, neighboring pixels tend to have similar values, you can use the cooccurrence matrix. The values in the matrix will be distributed mostly along the diagonal.

In the snapshot below, the second image is random and non-correlated. For that image the cooccurrence matrix is pretty much a constant.  I'm not showing the matrices themselves, but a plot of their values such that white means '0'.


derivation of discrete Fourier Transform

The Discrete Fourier Transform is one-and-the-same as the Discrete Fourier Series.  The DFT (and explicitly DFS) maps one periodic sequence, $x[n]$ of period $N$:

$$ x[n+N] = x[n] \qquad \forall n \in \mathbb{Z} $$

to another periodic sequence, $X[k]$, having the same period $N$:

$$ X[k+N] = X[k] \qquad \forall k \in \mathbb{Z} $$

and the inverse DFT (or "iDFT") maps it back.

I have a pretty rigid Nazi-like position about it (see the other answer), but I don't think that the math cares.  Now if you wanna connect this to the more common (but often misinterpreted) derivation, you can look at a previous answer of mine.

The term "support of $x[n]$" simply means the domain of values of $n$ in which $x[n]$ (presumed to be non-periodic and to have "finite support") is potentially non-zero.  It could also mean the duration of $x[n]$.  But when you periodically extend $x[n]$, then that sequence has no finite support.  $x[n]$ is defined and potentially non-zero everywhere (or all integers).


Lowpass : What is wrong with average filtering?

The question was : What is best between mooving average and 1st order filter when no RAM is available ?

The solution would be : None of them, except for a specific need : rejection of a particular frequency. Otherwise, it is better to use filter of higher orders than 1st order.

Mooving average :

Drawbacks : delay and poor rejection.
Assets : huge rejection at several particular frequencies

1-order filter

Drawbacks : rejection which is poor whatever frequency
Assets : no delay

Higher order filters. Higher order filters include : cascaded 1st order filters (here under an exemple with 5 cascaded 1st order filters), butterworth filters, etc.

Assets : better rejection than 1st order filter and no need of much RAM. Small delay.

View of the result :

Thank you @user883521 and @Florent Ecochard for your help !

Here is the Pyhton code for a 5th order filter (for other orders, cf links here under. for other filters, see signal.lfilter)

import numpy as np
import matplotlib.pyplot as plt

nfft = 8192
fsampl = 200
dtsampl = 1 / 200

fcutoff_1st_order = 1
dtcutoff_1storder = 1 / (2 * np.pi * fcutoff_1st_order)
fcutoff_5thorder =2.635
faverage2 = 3

print("generate signal")
raw_signal = np.zeros(nfft)
raw_signal[int(nfft/2)]=1

print("filter signal")
filtered_signal = 0.*np.copy(raw_signal)
for i in range(1, np.size(raw_signal)):
    filtered_signal[i] = filtered_signal[i-1]+ (raw_signal[i-1]-filtered_signal[i-1]) * dtsampl / dtcutoff_1storder
filtered_signal5= np.copy(raw_signal)
x = np.exp( -2 * np.pi * fcutoff_5thorder / fsampl )
for i in range(5, np.size(raw_signal)):
    filtered_signal5[i] = np.power( 1 - x, 5 ) * raw_signal[i - 1] \
                          + 5 * np.power( x, 1 ) * filtered_signal5[i - 1] + \
                          -10 * np.power( x, 2 ) * filtered_signal5[i - 2] + \
                          +10 * np.power( x, 3 ) * filtered_signal5[i - 3] + \
                          -5 * np.power( x, 4 ) * filtered_signal5[i - 4] + \
                          +np.power( x, 5 ) * filtered_signal5[i-5]
averaged_signal_66 = np.copy( raw_signal )
num_average = fsampl / faverage2
for i in range(1, np.size(raw_signal)):
    averaged_signal_66[i] = np.average( raw_signal[int( max( 0, i - num_average ) ):i] )

print("get FFT")
raw_signal_fft = np.fft.rfft(raw_signal)
filtered_signal_fft = np.fft.rfft(filtered_signal)
filtered_signal5_fft = np.fft.rfft(filtered_signal5)
averaged_signal66_fft = np.fft.rfft( averaged_signal_66 )
raw_signal_fftlog = 20*np.log10(np.abs(raw_signal_fft))
filtered_signal_fftlog = 20*np.log10(np.abs(filtered_signal_fft))
filtered_signal5_fftlog = 20*np.log10(np.abs(filtered_signal5_fft))
averaged_signal66_fftlog = 20 * np.log10( np.abs( averaged_signal66_fft ) )

print("plot")
freq_axe= np.linspace( 0, fsampl / 2, int( nfft / 2 ) + 1 )
plt.plot(freq_axe, raw_signal_fftlog, '-')
plt.plot(freq_axe, filtered_signal_fftlog, '-')
plt.plot(freq_axe, filtered_signal5_fftlog, '-')
plt.plot( freq_axe, averaged_signal66_fftlog, '-' )
plt.plot([0.02, 100], [-3, -3], 'b--')
plt.plot([2, 2],[-50, 0],  'k:')
plt.plot([3, 3],[-50, 0],  'k--')
plt.legend(['raw signal (dirac)', '1-order 1Hz filter', '5-order 1Hz filter', 'mooving average 66 samples',
            '-3dB', 'aliasing freq', 'unwanted signals frequency'],  loc="upper right")
plt.ylabel('amplitude (dB)')
plt.xlabel('frequency (Hz)')
plt.xlim([0,8])
plt.ylim([-30,2])

plt.show()

useful links :
http://www.analog.com/media/en/technical-documentation/dsp-book/dsp_book_Ch14.pdf
http://www.analog.com/media/en/technical-documentation/dsp-book/dsp_book_Ch15.pdf
http://www.analog.com/media/en/technical-documentation/dsp-book/dsp_book_Ch19.pdf


How to determine the frequency of a waveform?

Question 1: the sampling frequency would be $64$ Hertz
Question 2: the frequency indices range from $0$ to $63$

In the following code, I tried to mimic your data. At each run, a different observation is produced, based on a sine with a different amplitude distortion and offset.

Using the FFTR.m code that displays one half of the spectrum for real signals, and removing the average, you can get the following figure:

% Simulated data to mimic your signal
timeSamp = 1/64;
nSample = 64;
xAxis = linspace(0,1,nSample)';
dataClean = sin(2*pi*xAxis*6);
distAmpl = 1/10*( medfilt1(rand(nSample,1),3)+0.15);
dataDist = data.*distAmpl + 0.01*(rand+10);

% Half-axis FFT for a real signal
 [fftR,fftAxe] = FFTR(dataDist-mean(dataDist),timeSamp);

% Display
subplot(2,1,1)
plot(xAxis,dataDist,'.-');grid on;axis tight
xlabel('Time');ylabel('Amplitude');
subplot(2,1,2)
plot(fftAxe,fftR,'.-');
xlabel('Frequency');ylabel('Amplitude');
axis tight;;grid on


Yes it is.
Y(t), X(t) have the same periodicity periode for the sin signal.
The cos signal is not periodic.

That is of course assuming that u(t) is the unit step function with u(0) =1.


Can software distinguish between different musical instruments

What you want to do is timbre recognition.


Well, check out Wikipedia.

The ratio is defined as:

$$
F = \frac{\mathrm{between-group\ variability}}{\mathrm{within-group\ variability}}
$$

where, using your notation, they say that the between-group variability is:

$$
\sum_{i=1}^{M} (u_i - u)^2/(M -1)
$$

and the within-group variability is:

$$
\sum_{i=1}^{M} \sum_{j=1}^{N}( x_{ij} - u_i )^2 / (N - M)
$$


To get a panning effect I know that there is a need for a Circular Matrix. The angle is can be set as a default or I can let the user to select the angle.

How do you know this? Can you cite a source? Audio Panning is typically done with one input and two outputs. It's either done with constant energy or constant amplitude pan. Constant energy is similar to what the second column of your matrix would do if the angle goes from 0 to 90 degree with 45 degrees being "centered"


Creating an array whose DFT is all whole numbers

The DFT matrix has $e^\frac{-2 \pi j}{2^N}$ as coefficients. For $N >= 3$, these coefficients might have irrational numbers as real/imaginary parts (eg $\frac{\sqrt 2}{2}$ for $N=3$); so all these irrational terms will have to cancel out exactly in the DFT computation. There are ways of achieving that (repeating is one of them) ; but this keeps you away from the general case that your code has to test! For example, this will create symmetries in your data that might make an error in the bit reversal, or in the higher orders (3rd level or above) of butterflies undetected!

Another thing to take into account is that due to the pure floating point implementation, even if the expected result should be a whole number, your FFT code, especially on large sizes, might find a result like 0.99999999 (Heck, some whole numbers are not even representable as floats!) ; and this wouldn't be incorrect. Plausible on doubles, very likely of floats.

So you'd rather use very generic data (ie: that will yield some irrational numbers) ; and some unit tests that guarantees you that your results are within some tolerance of values computed by a reference implementation such as FFTW. For example in gtest there is EXPECT_FLOAT_EQ or EXPECT_NEAR to check for float equality (up to representation errors) and for float equality given a tolerance.


Signal Model

The first thing you really need when you start to use a Kalman filter is what I would call the "signal model".  That means, you need to take a guess (and sometimes that's all it is, though usually you can have an educated guess) at how the signal you're measuring, $y$, was generated.

The usual signal model comprises $F$, $G$, $H$, $Q$ and $R$:

$$
x_{k+1} = F x_k + G u_k + w_k\\
y_k = H x_k + v_k
$$

where $Q$ is the covariance of $w_k$ and $R$ is the covariance of $v_k$.

What are you trying to achieve?

Implicit in the signal model is that you've chose the "state variables" (the components of $x_k$) that you are working with.  Usually, what you are trying to achieve is to measure something and, from that, infer the value of one or more of the state variables.

What's missing?

Based on what you've said, the only thing that seems to be missing is how you get from the states $x_k$ to your measurements $y$ (i.e. $H$).


I think part of the problem is an awkward and inconsistent naming convention. There are 4 flavors of Fourier Transforms depending on which domain is continuous or discrete (which maps to being aperiodic or perodic in the other domain). So we have

         Name                  Time                  Frequency
Fourier Transform    continous/aperiodic     continous/aperiodic
Fourier Series       continous/periodic      discrete/aperiodic
Discrete Time FT     discrete/aperiodic      continous/periodic
DFT or DFS           discrete/periodic       discrete/periodic

A better naming would have been

         Name                  Time                  Frequency
Fourier Transform    continous/aperiodic     continous/aperiodic
Fourier Series       continous/periodic      discrete/aperiodic
Discrete FT          discrete/aperiodic      continous/periodic
Discrete FS          discrete/periodic       discrete/periodic

so that discrete refers to "discrete in time and periodic in frequency" and "series" refers "discrete in frequency and periodic in time". In other words "series" means sums and "transform" means integrals. Discrete mean sums and continuous means integrals.


Depth estimation from live video stream in matlab

Since you have two webcams you can reconstruct the video scene using the frame sequence. If your webcams happends to be stereo pair (with known extrinsic parameters aka rotation and translation) your task becomes pretty easy. On the other hand if webcam locations are unknown the problem gets a little trickier, but it can still be solved using matlab's computer vision toolbox.
In both of the above cases the scene reconstruction is accurate up to a scale factor. In order to calculate the depth in real units you need length knowlegde about an item in the scene. i.e the edge of a tracked square is 15cm.


npt is the number of frequency points that are used to define the desired frequency response. It's the length of the inverse FFT that is applied to the frequency domain data. It defaults to $512$ points, but if you want to design very long filters (with many taps) then you should choose a larger number.

lap defines the width of the transition band, if there is one. There is a transition band if there are two equal frequencies specified in the vector f, where a step in the desired response occurs. The wider the transition band the smaller the approximation error in the pass bands and stop bands, so there is a trade-off between steep transitions and approximation errors in the pass bands and stop bands.


Just take the real part of the complex signal.

(And if you do not believe it: click here.)


Complex Conjugate Sinusoids in Forward DFT

Really the DFT is just the inner product of your signal with complex exponentially of various frequencies.

Taking the inner product of complex vectors requires taking the complex conjugate.

Reading your code, I do not believe this is correct:

cp_sin = complex(generate_sine(k, 0, 1, Fs, N) / (N/2), generate_sine(k, pi/2, 1, Fs, N) / (N/2));

shouldn't it make a complex conjugate:

cp_sin = complex(generate_sine(k, 0, 1, Fs, N) / (N/2), -generate_sine(k, pi/2, 1, Fs, N) / (N/2));

where I've placed a "-" in front of the second generate_sine call.

Also, can you let us know what the arguments to generate_sine do. I'm not sure what the third and fourth arguments achieve?


Getting a weired behavior when windowing an audio signal

It turns out that the windowing function that is used performs by default "Zero-Phase Zero-Padding" on the input signal.
You can disable this feature by setting the zeroPhase parameter to False as follows:

w = std.windowing(type='hamming',
                  size=len(audio_pool['frames.actual-frames'][22]),
                  normalized=False,
                  zeroPhase=False)


Detecting trough widths and locations in 1d signal

A wavelet analysis or wavelet transform will give you scale(frequency) vs time vs amplitude. I am not sure how that will help you with your problem.
Here is an algorithm for wavelet transform :

Choose your wavelet : Haar , Molet etc ..say $ \psi (n) $
Choose Scale : the scale should be a multiple of smallest time interval dt.
Calculate the FFT of wavelet - $\Psi_0 (s\omega _k)$
Normalise FFT of the Wavelet - $ \Psi (s\omega_k)=(2\pi s / \delta t)\Psi_0 (s\omega _k)$
Find FFT of time series x(t) = $ X_k = 1/N  \sum_{k=0}^{N-1} x(n) e^{-2\pi i k n/N } $
Then Wavelet transform is :
  $ W(s)= \sum_{k=0}^{N-1} X_k\Psi^* (s\omega_k) e^ {i\omega_k n \delta t} $    ...where * is complex conjugate
Generally the Fourier transform of the wavelet should be known analytically.

Few good sources :

http://paos.colorado.edu/research/wavelets/bams_79_01_0061.pdf
http://users.rowan.edu/~polikar/WAVELETS/WTpart1.html


A clarinet has no even harmonics. What would produce no odd harmonics?

Producing "only" even harmonics is equivalent to producing all harmonics.

Say your original signal contains all harmonics:

1f, 2f, 3f, 4f, 5f, ...

Now remove all the odd harmonics:

2f, 4f, 6f, 8f, 10f, ...

The original 1f has been removed, so now you will hear 2f as the fundamental, and 4f will sound like the 2nd harmonic, 6f will sound like the 3rd harmonic, etc.  You're back to 1f, 2f, 3f, 4f, 5f, ... except an octave higher.  I think this might be something that happens in instruments, like "overblowing"?  But I'm not sure.  Anyway it will just produce a harmonic sound an octave higher, but with a different timbre.

Somewhat related: https://dsp.stackexchange.com/a/6178/29


About the order of a difference equation

To my understanding, the order of a difference equation is equal to the maximum delay, regardless if its delay on the input or output side. So, you have two first-order difference equations here.


Is interpolation of an audio signal to increase frequency resolution possible?

I do not understand why you have 30Hz resolution so I will focus only to the principle of the question "does interpolation increase resolution?".

Short answer is no, no new data, no new information.

A longer answer needs the spectrum visualization below, with time domain, continuous frequency domain and discrete frequency domain from left to right.

The interpolation technique here is to preserve the information of spectrum. DFT works on discrete frequency domain which is the part from $-f_s$ to $f_s$ of the continuous frequency version.

First, look at the continuous frequency domain, if you upsample your signal correctly, it is equivalent to changing the sampling frequency. You wish to double the number of data, but it is just removing one-half the spectrum replicas of sampling process.

Now, look at the discrete frequency one. This version is normalized from $-f_s$ to $f_s$ of the continuous frequency counterpart. $f_s$ is doubled, the spectrum is then shrinked by a factor 1/2. If we call $0 < \alpha < 1$ the proportion of non-zero frequencies before upsampling, this proportion is $\alpha/2$ after upsampling. Before upsampling, DFT gives you $N$ bins for $\alpha$ then $N\alpha$ bins for the spectrum; after upsampling it is $2N$ for $\alpha/2$ then always $2N \times \alpha/2 = N\alpha$ bins for the same spectrum. No, your resolution does not change at all.

To have smoother resolution, the only way is to add more data. In your example, instead of dividing to 30 frames, divide your audio file to 15 frames.


Periodicity and fundamental frequency detection algorithms

This is essentially the problem of Pitch Detection.  A fairly simple approach is to compute the autocorrelation (FFTs can be used to do this quickly),  and look for the peak with T > T_min