Fourier transform of function division in time domain

For division there is no equivalent to the duality between multiplication and  convolution. Note that from the existence of the Fourier transforms of $f_1(t)$ and $f_2(t)$, you cannot conclude anything about the existence of the Fourier transform of $f_3(t)=f_1(t)/f_2(t)$. So $f_3(t)$ might not even have a Fourier transform, and if it exists, it cannot be expressed in any simple way in terms of $F_1(\omega)$ and $F_2(\omega)$.


why does/ doesn't phase count when you analyze a voice signal?

Phase will have an effect on the mathematical analysis of speech signals however the traditional practice of speech (audio) processing has evolved based on the observation that human listeners do not have a phase sensitive hearing mechanism, which indicates (supported by experiments) that phase may be neglected in most speech processing applications, at least when compared to the frequency magnitudes. You may search the site for other related answers. But I think the most extensive source will be the audio engineering society.


Discrete Fourier transform of even function

The real question is: what did you expect to see and why? I give you a few ideas and hints to help you understand the problem:

if you use the time domain data as shown in the top figure, then the first (left-most) point is interpreted as time index $n=0$ (there goes the cosine ...). Just look at the formula for the DFT and you'll see why.
Second, the implicit periodic continuation of the time domain data in the DFT means that you don't have a pure sinusoid (just append the first data point after the last data point and so on, and you'll see the problem). What you probably wanted is an integer number of periods in your time window (= DFT length), because then the periodic continuation will not change your nice sinusoid.


The solution lies in the use of the x's. You need to convert the x's you use to lower frequencies, which can be either positive or negative. So only use x's between -N/2 and N/2 and the problem is solved.

def convertfreqs(x):
    return np.mod((x + N/2),N) - N/2

This works because it is the higher frequencies which cause the signal to change heavily inbetween pixels. Higher frequencies means higher x's.


Repeated Fourier transform - what happens?

if you define the continuous Fourier Transform in a unitary manner, my preferred unitary definition is

$$ X(f) \triangleq \mathscr{F}\{ x(t) \} = \int\limits_{-\infty}^{+\infty} x(t) \, e^{-j 2 \pi f t} \ dt $$

$$ x(t) \triangleq \mathscr{F}^{-1}\{ X(f) \} = \int\limits_{-\infty}^{+\infty} X(f) \, e^{+j 2 \pi f t} \ df $$

then you can see a lotta symmetry and isomorphy in the forward and inverse transformation.  in fact they are exchangeable since $-j$ and $+j$ both have equal claim to squaring to be $-1$.  (i.e. if, in all of our textbooks and technical and math lit, every $j$ was replaced by $-j$ and vise versa, all of our theorems would be just as valid.  the choice of which imaginary unit to go with which direction of Fourier transformation is a convention.)

now, it's not hard to see where the duality theorem comes from.  Given the above, then

$$ x(-f) = \mathscr{F}\{ X(t) \} = \int\limits_{-\infty}^{+\infty} X(t) \, e^{-j 2 \pi f t} \ dt $$

$$ X(-t) = \mathscr{F}^{-1}\{ x(f) \} = \int\limits_{-\infty}^{+\infty} x(f) \, e^{+j 2 \pi f t} \ df $$

from that, it's not hard to see that

$$ x(-t) = \mathscr{F} \Big\{ \mathscr{F}\{ x(t) \} \Big\} $$

and that

$$ x(t) = \mathscr{F} \Big\{ \mathscr{F}\{ x(-t) \} \Big\} $$

so it's not hard to see that if $x(t)$ is even symmetry ($x(-t) = x(t)$) then transforming it twice gets you back to the original.

and this

$$ x(t) = \mathscr{F} \bigg\{\mathscr{F} \Big\{\mathscr{F} \big\{ \mathscr{F}\{ x(t) \} \big\}\Big\}\bigg\} $$

so it's sorta like multiplying by $j$.  do it four times and you wind up with the thing that you started with.

you can use this fact to create an infinite number of Fourier transform pairs that are exactly equal to each other.  all you have to do is create the Fourier transform three levels deep and add each to the original.  each time you FT that, you get the same thing.

there are two simple functions that i can think of that have themselves as their own FT.  one is the Gaussian function and the other is the Dirac comb (both properly scaled).


You must remember that $s=j\omega$ and that $\omega=2\pi f$ and the transfer function is:

$H(s)={\prod_n (s-z_n)}/{\prod_m(s-p_m)}$

Where $z$ and $p$ are the zeros and poles respectively. Then you just need to adapt to your case.


Locating reference audio in longer audio

It seems that a matched filter is good enough for my use case.


Matlab: if given nothing but an audio file, how do I modify it to produce incrementally lower signal to noise ratios?

I like your noise model, I think that is really the best way to find SNR if you only have access to recordings. Using this same model couldn't you take your backgroundnoise scale it (duplicate it or trim it as necessary) and then add it to the signal? In this sense the noise wouldn't be entirely random but somewhat system specific

edit

in order to calculate your new noise value you would simply rearrange the SNR calculation solving for noise. Then scale the original noise accordingly. you can use some code like this

desired_snr_db = 10;    %change this value

est_signal = audio - backgroundnoise;

%the original matlab equation for SNR is:
%signalPow = rssq(x(:)).^2;
%noisePow = rssq(y(:)).^2;
%r = 10 * log10(signalPow / noisePow);
%solving for noisePow we get
noisePow = rssq(est_signal(:)).^2 / (10^(desired_snr_db / 10));

%assuming rssq(new_noise).^2 = (scale_factor * rssq(old)).^2
scale_factor = sqrt(noisePow / rssq(backgroundnoise).^2);

new_noise = scale_factor * backgroundnoise;

%add noise to our estimated signal
new_signal = est_signal + new_noise; %you may need to do stretching
                                     %trimming to get the right sizes

new_snr = snr(new_signal - new_noise, new_noise)

if you want you can double check the snr by finding a lull in the new signal and comparing it to the known sounds just as before


Percentage difference between two signals

At the least you should normalize it by the length of the sequence. But what you're trying to do is a pretty standard thing. I'm not sure what percentage difference really means. What you should probably be using is something like the root mean squared error:

$$RMSE=\sqrt{\frac{1}{N}\sum_{n=0}^{N-1}(x(n)-y(n))^2}$$

where $x(n)$ and $y(n)$ are your two sequences.


Simple analog anti aliasing filter good enough if signal frequency is far from sampling frequency ? What about equivalent time sampling?

The sampling theorem states that you need at least "two samples per one Hz of bandwidth". Anti-aliasing filters are required when the bandwidth of the original analog signal is not well defined and a finite bandwidth needs to be enforced.

In your application that doesn't seem to be the case. If all you want to sample is really just an unmodulated steady-state sine wave, the bandwidth or your original signal is already limited. No anti aliasing filter is needed and you can simply sample at any rate slightly larger than twice the sine wave frequency. 22.05 kHz or 24 kHz in your case would work fine.

ETS is a bit tricky since you need some sort of a trigger mechanism and then a variable delay off that trigger. The equivalent sample rate depends on your delay steps and the whole process depends on the signal being periodic with high precision.

You actually CAN sample at a much lower rate. The theorem requires 2 samples per Hz of bandwidth. For a pures sine wave the actual bandwidth is almost zero so you can sample with a very low sample rate and still don't have any loss of information. You just need to manage the mirror spectra properly, i.e. your application needs to be aware of a frequency offset.


How to discern the distribution function of mean-squared coherence values?

One approach might be an old method from "An approximation to the cumulative distribution function of the magnitude-squared coherence estimate" by Nutall and Carter (1981), who basically apply a non-linear distortion to transform the MSC to near-Gaussian random variables, then you could apply standard methods.

http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=1163657

A more recent approach used overlapping segments to create an (approximate) PDF and CDF for the MSC: "Approximation of statistical distribution of magnitude squared coherence estimated with segment overlapping" by Bortel and Sovka (2006):

http://www.sciencedirect.com/science/article/pii/S0165168406003471

They give a procedure for the comparison of two MSC estimates.

There may be more refined methods, but I'm not aware of them.


Forward/reverse median filter. Reference needed

I've held off on answering this because I don't have MATLAB to test any of this.

I've employed the same trick with exponential smoothing.  The purpose of this is to have the lags in the opposite direction cancel each other out, so the resultant smoothed signal is centered instead of shifted from the source signal.  For sinusoids, the lag cancellation is exact.

So, I thought this might be the same motivation in this.  It does not seem to be the case.  (Side note:  In my searching I found the MATLAB function "movmedian".  I don't immediately see how it is that different from "medfilt1".)

I agree with Laurent Duval that if an odd number is used that the forward and backward filtered signals will be the same.  However, OP is using an even number.  It seems to me that the backward filtered signal will be just one index off of the forward one since the window size is the same and the median for each window should give the same result.  Thus, excepting the end cases, you should be able to just do the forward one and then define:

$$ M[x] = ( Mfw[x] + Mfw[x+1] ) / 2 $$

Which is a very simple FIR.

Two cents minus inflation worth.

Ced


Does lowpass filtering and decimating improve SNR?

So, let's not forget what SNR is: it's a relation of powers present.

The thing that improves SNR is a propoer low-pass: it leaves the signal power alone and reduces the power of the noise.

An ideal low-pass filter will leave zero noise outside its specified bandwidth – so it doesn't matter whether you "cut off" these bandwidths using decimation or not, from a pure SNR perspective.

  This question is of practical importance: when measuring a phenomenon, say at 1Hz, is it useful to use a very high sampling frequency such as 1kHz or more in order to improve the SNR by filtering and decimating (given that the memory and computational cost is not a problem, of course)?

Yes! That is called oversampling, and it's useful for a couple of reasons, foremost for doing so-called noise shaping, which allows you to get a lower noise floor in your band of interest. With appropriate filtering, you do get a lot of SNR gain.

This is actually done pretty often; Software Defined Radio devices often run at dozens to hundreds of Megasamples, no matter how much bandwidth you want to observe, and do the filtering and decimation in-device. That way, you get the oversampling benefits, and also, you don't waste money and potential phase linearity etc by implementing steep-edged analog anti-aliasing filters, but can get away with much more relaxed filter specs.


Symbol synchronization with unknown symbol rate

The one idea I have is to first use a symbol rate estimator, but the problem is that there will be some error from the estimator, then the error gets used in the symbol synchronization. Ideally, there is some synchronization method that don't need to know the symbol rate ahead of time or one that is resistant to small errors in the symbol rate estimate.

I mean, you need some tolerance anyways. No two oscillators in this universe are exactly the same frequency. And oscillator price pretty drastically scales with required long-term accuracy... so, you must assume that your receiver doesn't have perfect symbol rate knowledge, anyway.
Symbol timing recovery methods are typically resistant to small rate errors – look at the early/late gate, especially: You'll just see that you're a bit early more often than late if you're overestimating your symbol rate.
You just need to continuously run that.
This really isn't much different to phase recovery methods, which typically can also correct small frequency errors.

How could a multi-symbol rate system use an early late gate?

Uff, in such a system you'll often use a bit more complex synchronization methods, which take information from preambles, and, just as mentioned above, will keep tracking the symbol rate.


Since your input is discrete the spectrum is periodic with the DFT length. $X_0$ is your DC term, but so is $X_{512}$, $X_{1024}$, $X_{-512}$, etc.
It really doesn't matter which exact period you pick: [0 511], [-255 256], [-256 255], etc. Any choice will work and you always get the same amount of data. Note that $X_{-256} = X_{256}$.


Does converting an FIR filter from causal to non-causal remove phase delay?

Correct. The magnitude response would be the same and the phase response would be summed with $\omega$ (omega).


$D$ is $m \times m$, where $m$ is the total number of time frames in the spectrogram, it represents the squared euclidean distance of each frame to each other frame. In Matlab you could calculate it like this:

% Equation 10
for k = 1:m
    for l = 1:m
      D(k,l) = sum((X(:,k)-X(:,l)).^2);
    end
end

This is clearer from the text than the equation itself.

The $P$ matrix is calculated for each frame from the $p$ the nearest frames to it. The $Y$ matrix is then the median of $P$ at each frame. In matlab you could calculate it like this.

[~,I] = sort(D,2);
for k = 1:m
   pNearest = I(k,1:p);
   P = X(:,pNearest);
   Y(k,:) = median(P,2);
end


Why can Sum of Squared Differences (SSD) seen as cross-correlation?

Short Answer:
The statement "the the squared terms are constant" is true only if some further theoretical assumptions are made. In practice, the terms are often not constant (i.e. they are dependent on $m$ and $n$). For many discussions however, it's ok to treat the terms as constant.
Details:
The term $\sum_{i,j}\left[g(i, j)\right]^2$ is not interesting since it is clearly independent on $m$ and $n$. The term $\sum_{i,j} \left[t(i-m, j-n)\right]^2$, however, formally depends on $m$ and $n$, as $m$ and $n$ appear in the summed expression. Nevertheless, it may be effectively constant depending on how the bounds for the running variables $i$ and $j$ are chosen and whether $t$'s support is bounded or not.
In theoretical discussions, the image function $t(x,y)$ is usually assumed to have bounded support. This means that there are minimum and maximum values of both $x$ and $y$, outside of which all image pixels can be assumed to be  zero. This assumption is plausible because a real digital image only has a finite number of pixels, e.g., $1024$x$768$ pixels.
The sum is often assumed to cover the entire space, that is $i$ and $j$ go from negative to positive infinity: $\sum_{i,j}$ really means $\sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}$.
Only if these two assumptions (bounded image and infinite sum) are silently made, then the squared terms are independent of $n$ and $m$. No matter what $n$ and $m$ are you will always "hit" all the pixels of image $t$:
$\sum_{i,j} \left[t(i-m, j-n)\right]^2 = \sum_{i,j} \left[t(i, j)\right]^2, \textrm{for all } n, m$.
In practical template matching, the sum will often cover only the pixels of the template, for performance reasons. This means that, in the (larger) image, only the pixels are regarded which are currently covered by the template given the current $m$ and $n$. In these cases, the squared terms are not independent if $n$ and $m$.
As a last note, in your squared difference implementation you can choose to omit calculating the squared terms and only regard the mixed term. This way, you get the same result as if you would have taken the infinite sum (up to a constant), and thus the theoretical considerations will hold again.


What subclass of nonlinear systems can be represented by Volterra series?

The system must be time invariant and smooth in the functional derivative sense. That doesn't guarantee that the Volterra series converges (like with Taylor series, there are pathological counter examples), but almost all systems that have these properties have a convergent Volterra series.

The problem in practice is that the required expansion order for good approximations of your system is often far too great to provide for a useful implementation of the system. A condition for quick convergence is, that the functional derivatives are well bounded and decay quickly with order.

A little more detail:

Consider a non-linear system $f:V\to W$, where $W$ is an inner product space with induced norm and $V$ is a normed vector space. With this, look at a map

$$\begin{align}
W\times V&\to K_W \\
(w,x) &\mapsto \langle w,f(v+\lambda x)\rangle
\end{align}$$

This map is anti-linear in $w$ and non-linear in $x$. We treat $\lambda \in K_V$ and $v \in V$ as parameters. Further below, $v$ will be the center of the Volterra expansion and $\lambda$ will reveal the different expansion orders.

If we hold $w$ constant, we get a functional

$$
\begin{align}
F_w:V&\to K_W\\
x&\mapsto F_w[x]=\langle w,f(v+\lambda x)\rangle
\end{align}
$$

For bounded $x$, this functional can be expanded into a series of powers of $\lambda$ around $v$:
$$F[x]=\sum_{n=0}^\infty \frac{\lambda^n}{n!}F^{(n)}(x)$$

where $F^{(n)}(x)$ is the $n$-th (functional) derivative of $F$ around $v$ in direction $x$.

The assumption that all $F^{(n)}$ exist and the series converges for all finite norm choices of $w$ is what we mean, when we say: "the system $f$ is smooth".

Going from here to the full Volterra series is not too difficult. The key is understanding, how the functional derivatives relate to $w$ and $x$. Since $F$ is anti-linear in $w$, every single term in the series expansion must be, and therefore $F^{(n)}$ as the only factor with dependence on the choice of $w$. It is also clear by the option to absorb $\lambda$ into $x$, that the derivative $F^{(n)}$ must be of $n$-th order in $x$.

Hence we can write $F^{(n)}(x)=\langle w, G_n(x) \rangle$ with multilinear $\hat{G}_n:V^n \to W$ as well as an embedding $E_n:V \to V^n$ with $E_n(v)=(v,v,\dots,v)$ composed to $G_n=\hat{G}_n \circ E_n$.

In infinite dimensional vector spaces with uncountable bases, these multilinear maps are expressible through nested integrals over kernel functions. If you use the time basis for this expansion and make use time-invariance, you get nested convolution integrals and arrive at the Volterra series. I'll leave this mostly notational part to you. If you need help with it, I can also elaborate more.

As the last step, you use the properties of the inner product to draw $w$ out of the series expression and simply drop it, as the resulting equation needs to hold for all proper $w$.


In the $2.4 \textrm{ GHz}$ band, the data rate is indeed just $250\textrm{ kbit/s}$ but the direct-sequence spread spectrum (DSSS) coding results in a bandwidth of $2\textrm{ MHz}$. This means that there will be flat fading only if the channel's coherence bandwidth is significantly greater than $2\textrm{ MHz}$, say around $10\textrm{ MHz}$ (this is just a fuzzy estimate). This corresponds to a maximum delay spread of approximately $1/(10\textrm{ MHz})=100\textrm{ ns}$.

Consequently, any environment with a delay spread greater than approximately $100\textrm{ ns}$ will result in frequency-selective fading. So in relatively small indoor areas the fading will be mostly flat, but in large office spaces or outdoors you can expect frequency-selective fading.


Image Processing and applicability of 2D Fourier Transform

Please help me understand why in image processing, transformation along both axes is needed.

A one dimensional signal describes how does a quantity vary across, usually, time. Time, commonly represented by the symbol $t$, is the only parameter required to describe completely the signal at $t$.

A two dimensional signal describes how does a quantity vary across two parameters that are absolutely required to describe it completely. When referring to images, the quantity that is described is, usually, radiant flux. That is, how much "light" (more generally, radiation) is received by the sensor. In common handheld cameras, each pixel of an image describes how much visible light is received by some point in the scene that is viewed.

The complete set of pixels of an image, describes the variation of visible light across the surface of the camera's sensor.

When applying the Fourier Transform to a one dimensional signal, the dimension of time is transformed to a dimension of frequency and the transform breaks the signal down to a sum of sinusoids.

When applying the Fourier Transform to a two dimensional signal, its two spatial dimensions are decomposed into sums of orthogonal spatial sinusoids. To cut a long story short, if the basis functions of the Fourier Transform were not orthogonal, this trick of decomposing and recomposing would not be possible. How does this look like? It looks like a carton of eggs:

The higher the spatial frequency, the smaller the eggs (more of them fit in the same length) and vice versa.

More formally, it looks like this:

And it doesn't even have to be "symmetrical", that is, each one of its dimensions may be supporting a different spatial frequency:

In this last image, we have many more sinusoidal cycles across the $x$ dimension than across the $y$ dimension.

Therefore, whereas in the one dimensional case, a signal as complicated as the voice of a singer is decomposed into a set of simple "whistles", in two dimensions, an image as complicated as Lena is decomposed into a set of elementary little blobs. In the first case, the signal is correlated with a series of progressively increasing frequency sinusoids, in the second case, exactly the same thing happens only now, the signal is a patch of pixels and the "sinusoid" is a patch of spatial frequencies that could vary differently across the $x$ and $y$ dimension.

Now, in terms of expressing this process via the one dimensional Fourier Transform (that performs this correlation process with one set of sinusoids), the same, is applied twice.

Consider Lena. And apply the Fourier Transform across each one of its rows. What do you get? You get a set of rows in the frequency domain. These describe the set of sinusoids that visible light across the rows of the image varies. But!!! at this point, we know nothing about the set of sinusoids that describe visible light variation along the vertical.

Another way to "visualise" this, is to consider the DC bin of the row Fourier Transforms (frequency is zero). This tells you the average brightness of the pixels in each row but it still varies along the column direction! (i.e. we know nothing about the DC along the columns).

In other words, where we had an $f(x,y)$, we pass it through the Fourier Transform along the rows and we now have an $f(F_x, y)$. We are now in an intermediate state where one dimension is frequency and the other is still space.

For this reason, you apply the Fourier Transform once more along the columns of the Fourier Transform of the rows. (And to get back to the DC example, you now get one DC coefficient that describes the average brightness along the rows and the columns, that is, you get the average brightness of the image.)

Now, remember, the one dimensional Fourier Transform decomposes a signal ($x(t)$) into two series of "strength" coefficients, one for the strengths of the $\sin$ and one for the strength of the $\cos$ at different frequencies. The two dimensional Fourier Transform does exactly the same thing but now the "strength" coefficients are two dimensional. That is, some coefficient at pixel $i,j$ (in the Fourier Transformed image, i.e. the Frequency Domain) describes the contribution of a "carton of eggs" with different number of cycles in the $x$ and $y$ dimension to the whole image.

Generalising to three and higher dimensions is done similarly.

Hope this helps.

(Note, all images retrieved via Google Images and linked, rather than uploaded, to the post)


What features to use for multilayer perceptron model of an EEG signal?

Fourier analysis doesn't provide you with temporal information.  If you need temporal information you have several choices, two of which are the short term Fourier transform STFT, and wavelet transforms.

The STFT provides frequency coefficients in fixed time intervals across your data.  There is a trade off between the time and frequency resolution obtained.

Wavelet transforms pick apart features in your data at various time resolutions which may give you better time localization for your coefficients.  Additionally, you can choose different wavelet functions to control how features of your signal are extracted.  Another way of saying this is that your choice of wavelet function affects what your coefficients represent.

Most of the work on feature extraction that I've seen recently (this is in general, no specific to EEG signals) employ some form of wavelet analysis.

Here is an example paper: http://uprt.vscht.cz/vav/RSoftware/S7help.pdf


Most accelerometers suffer from bias --- a non-zero, constant rest measurement. See, for example, this page.

  A further imperfection is called zero g bias level, or offset. This is the reading the sensor provides with no acceleration, and, unfortunately, that is rarely zero (in the Z direction it should be 1 g). In fact, the error can be substantial, frequently hundreds of times larger than the device’s resolution.

So you need to remove this DC offset (which is how the bias appears) from all your measurements.  If you need to know the $g$ value from the $z$ accelerometer, you will need to calibrate it first in the $x$ or $y$ plane (to remove the effect of gravity).

Have a read through this answer to another question, and try applying a DC blocker to your data.


Is it possible to observe more area in azimuth than physically traveled in stripmap SAR?

It's so that the cross-range (along-track) resolution is consistent for all targets in the scene. There is no reason why you could not use the flight path between [-5,5] to form the same size image. However, in that case, the targets on the edge of the scene will have been illuminated for a smaller angle than the targets in the middle of the scene, which will given them a smaller region of support in the frequency domain and, hence, a coarser resolution.

When choosing the first and last sensor locations that are used to form the image, it's common to select the first point to be that where the beam first enters the imaged scene and the last point to be when it leaves it. Picture the footprint of the beam on the ground as the sensor flies by; it's usually an oval or cone or something like that. At -20 is where the edge of the oval (or whatever shape it is) begins to intersect the imaged scene. At +20 is where it stops intersecting it.

I should add that, in general, the size of the imaged scene is limited only by the PRF of the sensor.


making pink noise (1/f) using list of frequencies

Pink noise, as you know, has contributing amplitudes inversely proportional to each corresponding frequency. It appears more sinusoidal than white noise (constant power spectrum) because the low frequency components have large amplitudes. The high frequency components do not contribute that much to the overall noise amplitude.

To recreate pink noise in MATLAB, I would use the following technique:

1.) Generate a square wave (ones(n) will create a square wave of width n).

2.) Take the Fourier Transform of the signal (fft(x)). This will give you the sinc function. The DC component will be 1.

3.) In a for loop, add noise to the amplitude of each frequency of the FFT sinc function. The range of values of the noise will be dependent on 1/f.

4.) Take the inverse Fourier Transform of the noisy FFT. Then subtract 1 to yield a pink noise vector of length n centered around $0$.


Strange behavior from frequency to angular frequency in FFT

This is probably the single most asked question on this forum. A similar one was asked  and answered just a few hours ago: Spectral leakage from mathematical point of view
Your expectations are wrong: you only get a spectral dirac impulse if the frequency of the complex exponential on he FFT grid, i.e. an integer multiple of the sample rate divided by the FFT length. In your case you have $N_{FFT} = 10000$ and $F_s = 1000$, so your FFT bin spacing is 0.1 Hz. 50 Hz is an integer multiple of 0.1Hz but $50/(2\pi)$ isn't and so the energy will spread out over the adjacent bins.
A different way to look at it: the DFT assumes that the time signal is periodic with the FFT length. Since your exponential does NOT have an integer number of periods inside one FFT frame, the periodic repetition creates a discontinuity at the frame boundary


Physical significance of resolution in image processsing

So the question really is, "what is a pixel?"

A pixel when you access the data/sample of the image is just a number. The capturing device has generated this data represents the visible area on a particular surface and the color of the pixel represents the light that is falling on that surface.

Hence, when resolution of the capturing devices is high, you are able to capture equivalent of finder details from the incoming source of light and its' variations. when the same image is brought down to a resolution of say 1/4 which means the resultant new image represents the average color/intensity of that 4x4 window around that point. Hence, the lower resolution more "pixelating effect" do you see.

There is one more aspect to see here. The details in image also represents the size of the object subject to the scene geometry. So for example, if you capture a shot in camera which has one flower and a leaf, we identify that if the leaf is roughly taking twice as more pixels than flower, then leaf in real life is twice as bigger (from that perspective) than flower. However, when zoomed out, and pointing to the macro scene, the same number of pixels could be occupied by a full tree. So we can find the relative dimension of the object based on their area in the image because image is nothing but a projection of that scene on the canvas. It is possible to know, however, the exact dimension of the object as well if we can apply some known distances on the canvas or if we know the distance of the object from the camera film.

So in essence, pixel represent the light intensity or the color of the object at a fixed grid on the scene - the finer the resolution - more pixel reprsent the same object in finer details and more information. And, distance between the pixels represents the relative distance between two points on the scene.


Beginner's question on LPF, Need for Sinc and windowing

Interesting question. So you are asking why we do not use FFT instead of commonly used filters (in your case Since function).
In terms of computational complexity, FFT wins over digital filtering. However though, zeroing out FFT bins, instead of filtering is not a good idea. There are many reason for that and high among them is the artifacts that are produced by FFT and FFT leakage.
FFT bins in other perspective can be seen as representation of  discrete sines and cosines. Your real frequency resolution (number of bins) depends on your signal length, the shorter your signal is the fewer bins are there to represent sines and cosines that embed in your signal (you might consider taking FFT of too long windows? the problem is then you must assume your signal is stationary). On the other hand, if your signal's sinusoidal components have frequencies other that frequencies represented by FFT bins (so called integer frequencies), you'll frequency bins will smear over their neighbors as below:

So when you zero out a bin like 400th bin and its neighbors in above image. You cannot remove 400th bin completely because some portion of it lives in other far located bins and also some energy of other bins (like 500th bin) smeared over here (so removing 400th bin will effect the component in 500th bin).

The other reason is that, you usually do not have your entire signal recorded and stored in memory to perform FFT, your data usually is streaming, and in this case only common methods of filtering come in handy, because they perform filtering on the data on the fly, whereas FFT requires all of the samples to be available. One might ask, what if we segment our data to small windows and then use, FFT? In that case, the leakage problem becomes severe and your filters accuracy degrades.

More interesting post on this regards can be find here(s):

https://electronics.stackexchange.com/questions/109820/why-should-i-use-digital-filters-for-bandpassing-rather-than-simply-manipulate-s
Why is it a bad idea to filter by zeroing out FFT bins?


Hi: I don't know whether white noise can have a non-zero mean ( Dilipe can tell us that )  but, if it can't, then just call it whatever you call white noise with a non-zero mean.
But, assuming it's zero, this is how you can do it. ( If it's not zero, then the  solution is similar but the last step would be a little different. Use
$$\operatorname{Var}\left\{X\right\} = E\left\{X^2\right\} - E\left\{X\right\}^2$$
Note also I wrote it with "t's" and "epsilon's"  just because I'm more used to those but you can think of $t$ as $n$ and $\epsilon$ as $w$.
\begin{align}
x_t &= 0.2 x_{t-1} + \epsilon_{t} + \epsilon_{t-1}\\
\implies x_{t}(1 - 0.2L) &= \epsilon_{t} + \epsilon_{t-1} \tag{$\scriptsize{\text{using lag operator}}$}\\
\implies x_{t} &=  \frac{\epsilon_{t}}{1 - 0.2L} + \frac{\epsilon_{t-1}}{1 - 0.2L}\\
\implies x_{t} &= \sum_{i = 0}^{\infty} 0.2^{i} \epsilon_{t-i} + \sum_{i = 0}^{\infty} 0.2^{i} \epsilon_{t-1-i}
\end{align}
So, since the mean is zero, the $E\left\{x^2_{t}\right\}$ is the variance of the RHS.
So, it is
\begin{align}
\sum_{i = 0}^{\infty} \left(0.2^{i}\right)^2 \cdot \sigma^2_\epsilon + \sum_{i = 0}^{\infty} \left(0.2^{i}\right)^2 \cdot \sigma^2_\epsilon&= \sum_{i = 0}^{\infty} 0.04^{i} \cdot \sigma^2_\epsilon + \sum_{i = 0}^{\infty} 0.04^{i} \cdot \sigma^2_\epsilon\\
 &=  \frac{2 \cdot \sigma^2_{\epsilon}}{(1- 0.04)}
\end{align}


Note that $\delta(t)$ is no function but a distribution. This means that it doesn't have function values but it is only defined by its integral. Consequently, the (generalized) function $x(t)=-\delta(t)+u(t)$ has no value at $t=0$ (the only useful value would be $-\infty$ but that's kind of hard to sketch). For sketching the function graph, you can define function values for $t\neq 0$ in the usual way. At $t=0$ you simply sketch a negative $\delta$-impulse, for $t<0$ the function is zero, and for $t>0$ you assign it the value $1$.


why does the decision boundary change in BPSK? is it true?

Zero-mean noise by itself can't modify the decision boundary. However, the number of things that can go wrong in your system is large.

Is your channel flat or frequency-selective?

You can think of the channel as the sound card and its driver, the loudspeaker, the air, the microphone, and the receiver's sound card and its driver. In my experience, there's very large variability in the frequency response of such a system. You can find the channel's response (for example by transmitting noise or an impulse) and use only the frequency band where it's flat. Or, you can design a wide-band system, but then you'll need an equalizer in the receiver.

Is your carrier synchronizer working?

Plot the received signal and see if its envelope is constant. If it's changing over time, you may have a faulty synchronizer. One quick way to get around this problem is to use DSB-LC modulation and demodulate the signal with an envelope detector.

Are you using a matched filter?

You need one to get optimal signal-to-noise ratio. After your filter, plot the signal's eye diagram. If it's open, then you know that your system up to this point is correct.

Is your symbol timing correct?

Since the clocks of the transmitter and receiver are different, you need to do adaptive symbol timing recovery to properly sample the matched filter's output. To test this, plot the sampled symbols over a long time. If they are correct for some time, but then start to get closer and closer, your symbol timing recovery is off.

Is your frame synchronization correct?

Finally, you need to verify your frame synchronization is correct; otherwise it's very hard to align the transmitted and received streams for bit error calculations.

I'd suggest verifying the correct operation of your system block by block, in the order I have listed them.


Identify the frequency spectrum of an audio file by using matlab

It looks like you are trying to recreate the functionality in spectrogram. Spectrogram will give you a short-time Fourier transform using whatever blocksize you specify.

Much of your code looks like it is re-inventing this.


Fourier transform is an isomorphism...but we don’t get when each frequency appears?

It's true that taking the Fourier transform will leave you without any (visible) information on time and vice versa, but of course you don't lose any information, you just represent it in a way such that in one domain you only see time information, and in the other you only see frequency information.

Take as an example the Fourier transform of a time-inverted (real-valued) function $x(-t)$: its Fourier transform has the same magnitude as the Fourier transform of the original function $x(t)$. The difference between these two Fourier transforms lies exclusively in the phase. So you're right to assume that timing information is encoded in the phase of the Fourier transform.

There is no time localization in the Fourier transform because its basis functions are complex exponentials extending from $-\infty$ to $\infty$. There are other transforms that will give you a certain degree of time and frequency localization, the most well-known of which is probably the Short-time Fourier transform. Also take a look at this related question and its answers.


Polyphonic detection/Mulit-pitch detection/Chord recognition

As @jan pointed out, you're probably asking for a little too much, especially if you're looking for a ready implementation. Doing a quick search on Google, I came across several papers that may be a helpful start.

In this paper called Multi-pitch Detection Algorithm Using Constrained Gaussian Mixture, the authors use the Expectation Maximization algorithm to solve a Gaussian Mixture Model to detect multiple tones. This is most likely not a very computationally efficient algorithm, since mixture models are normally very hard to solve, and have to be primarily done offline (no real-time version).

If you're looking for a moderately robust, but potentially fast algorithm, look at this Stack Overflow answer. The author recommends MUSIC and ESPRIT algorithms.

The best thing to do is to give something a try and come back to us with specific DSP and algorithms questions and sites like Stack Overflow for specific implementation questions.


How do I play and record audio at the same time?

I suspect that your problem is to do with having your speakers and microphone on different ALSA audio devices or subdevices.

You need to find the available sound cards for recording and playback on your system. The following commands will tell you this information :

aplay -l
arecord -l

You will need to define a default full duplex sound card using your choice of devices you see in the output to the above commands.

The way to define this full duplex card (made of two different soundcard devices) is to define it in your ~/.asoundrc file. This example shows you how to do that, however I will repeat it here.

The requirement is that you use plughw as your device name because the clocks on the two sound cards will drift and you will need resampling of some form. But you can play with this option.

So for example if you want to take input from a device on the first card (hw:0) and output to the second device on the second card (hw:1,1) then you would write the following in your ~/.asoundrc file :

pcm.!default {
    type asym
    playback.pcm "plughw:1,1"
    capture.pcm  "plughw:0"
}
ctl.!default {
        type hw
        card 1
}

Hopefully now Matlab sees a full duplex card as the default device and allows full duplex operation.


Why is it easier to extract color in HSV model than in RGB?

Say you have an object that is blue, and you adjust the brightness of a lamp in the room, or you move the object so that it does not get as much light. In RGB, all the three components (red, green, blue) will change but in HSV mostly the V (value) component changes. So to recognize pixels belonging to the object by its color, you only need the H (hue) and S(saturation) components of HSV, in contrast to needing all three components of RGB.

Well actually for almost black pixels the H and S components can be anything due to noise, so you will need to see that V is high enough. An equivalent check needs to be done with RGB too.


Differentiating between zero and not sending for OOK

A simple way to do it would be using a start bit like UART. When idle, don't transmit anything. When you want to transmit something, first transmit a start bit which is a '1' (ON). Then transmit a fixed number of bits known to both the transmitter and receiver. The receiver looks for the start bit and knows to begin demodulating after it receives the start bit.

A more rigorous way to do it would be the following:

Assemble bits into a frame of size X bits.
At the beginning of the frame, transmit a preamble sequence. This is some known sequence of bits that will be the same every time. This is what is done in Wi-Fi, Ethernet, etc.
The receiver looks for that preamble sequence by correlating the received signal with the expected preamble signal. Once you get a high peak in the correlation, the receiver knows that is the start of the frame.
The receiver is now synchronized and knows that it needs to demodulate X bits. After that, the receiver goes back to looking for the start of the next frame.

Raw modulation isn't really useful by itself; you need to add some extra stuff to have a practical communication system.


fftshift in MATLAB with even number of data points in double sided spectrum

If the size $N$ of the DFT is even, only one "extremal" point (after fftshift) is Nyquist.  If $N$ is even, you cannot have an arbitary $c_{-5}$ and $c_5$.  They must add to be whatever your $c_{-5}$ term is.

If the input to the DFT is purely real and if $N$ is even, consider the Nyquist point, at $c_{N/2}$ to be split in half.  One half is the negative-frequency component and the other half is the positive-frequency component.  Since the negative and positive frequency components must be complex conjugates of each other, and since $c_{-N/2}=c_{N/2}$ after splitting, that means that $c_{N/2}$ must be purely real just as $c_0$ must also be real.


Huffman coding for $8\times 8$ blocks of an image

Huffman coding, as an example of statistical entropy coding, is used as the last step of still image compression to further reduce any statistical redundancy (possibly) present at the output representation of quantized DCT coeffcients of a given block of N by N pixels (usually N = 8)

When the statistical behaviour of those DCT coefficients are analysed throughout the whole image, it will be seen (for naturally generated images) that certain symbols occur more frequently to the extend that using a variable length code (VLC) will become benefical instead of a fixed length one.

Thus to utilize this further possibility of statistical compression, VLC is employed via Huffman code tables for each of the AC, DC, Gray and Color channel coefficients separately according to their respective statistical symbol frequency distributions.

Now, the proper way of generating a statistically optimized Huffman table per block, actually necessitates a per-block optimized huffman table as you outlined in your question. However due to several reasons, one of which is again mentioned as the complaint in your question, such a method would oppose the initial purpose of image compression; to reduce its size.

Mainly:
1- inlcusion of per block optimized huffman tables into each block will reduce benefits of compression (it can even increase the block size for high quality compression where there will be so many symbolls and the table will include as many data as the block itself)

2- The loss of the benefit from discarding a per-block optimized Huffman table utilization will be mostly insignificant (i.e. A single whole image optimized Huffman table will perform almost as good as a per-block optimized one - even when excluding the effect of the table storage per block)

3- It will be more robust to create the table from an ensamble of blocks rather than a single one. (However note that due to statitical nonstationarity of typical images, the statistical behaviour from one block to an adjacent one other will change and hence the ensamble is not always a qualified source of data for table optimization)

4- Nevertheless, the approach taken is to use a single Huffman table for the whole image and store this inside the image file headers.

5- Now this single Huffman table can be defined in one of two ways: as either a pre-defined standard table or as an image specific optimized table. Most often the former approach is taken as the benefit from using specific image optimized Huffman table is attenuated by the table processing time and the complexity of programming required for per-frame Huffman table optimization.


What does the variable $\tau$ mean in frequency modulation?

It's just a dummy integration variable. The time variable $t$ is the upper integration limit, so you can't use it as the integration variable. An expression like
$\int_0^tm(t)dt\qquad(???)\tag{1}$
doesn't make sense, whereas
$$\int_0^tm(x)dx\tag{2}$$
does.
Take a look at this question and its answers over at math.SE.


It probably has the potential to do so, but the probability is marginal. Result of equalization depends on the initial pixel intensity distribution.

Image histogram equalization can have different implementations.
In one possible implementation, you would re-map a whole block of a given intensity level into a new calculated intensity. In such a case, objects will not change shape; boundaries will remain the same.

The aim of equalization is not necessarily the creation of a new image with strictly flat histogram will all bins equally occupied. That won't create a good looking image in any sense. Instead, the aim is to expand the pixel intensities to fully utilize the available dynamic range in the hope of getting a higher constrast with more visible dark regions and unsaturated highlights.


derive an average waveform from an FFT spectrum, where all cycles have same shape and period

I found the answer here:

https://www.dsprelated.com/thread/3253/simple-and-accurate-phase-estimation-of-harmonics-in-a-signal

Just change atan() to atan2() to account for the correct phase in all quadrants.


Let's have a look on the following model:

$$ y \left[ n \right] = \left( h \ast x \right) \left[ n \right] + \left( g \ast w \right) \left[ n \right] $$

Where $ x \left[ n \right] $ is the signal of interest and $ w \left[ n \right] $ is the AWGN with unit Variance.

In Matrix form it is written by:

$$ \boldsymbol{y} = H \boldsymbol{x} + G \boldsymbol{w} $$

Where $ H $ and $ G $ are the convolution matrices of the model.

Let's define $ v = G \boldsymbol{w} $ then $ v \sim \mathcal{N} \left( \boldsymbol{0}, G {G}^{T} \right) $ which implies (If the input to Linear Operator is Gaussian Variable then the output is also Gaussian Variable):

$$ \boldsymbol{y} \sim \mathcal{N} \left( H \boldsymbol{x}, G {G}^{T} \right) $$

Then the Maximum Likelihood is given by:

$$ \arg \max_{ \boldsymbol{x} } \det \left( 2 \pi G {G}^{T} \right)^{-\frac{1}{2}} {e}^{ -\frac{1}{2} {\left( H \boldsymbol{x} - \boldsymbol{y} \right)}^{T} {\left( G {G}^{T} \right)}^{-1} \left( H \boldsymbol{x} - \boldsymbol{y} \right) } $$

It is easy to see you only need to deal with the term in the power of the Exponent, so it is equivalent of:

$$ \arg \min_{ \boldsymbol{x} } {\left( H \boldsymbol{x} - \boldsymbol{y} \right)}^{T} {\left( G {G}^{T} \right)}^{-1} \left( H \boldsymbol{x} - \boldsymbol{y} \right) $$

Which is nothing more than Weighted Least Squares problem which matches the result of Linear Regression for Colored Noise (As expected for Gaussian Noise).

The optimal solution is given by:

$$ \hat{x} = {\left( {A}^{T} {\left( G {G}^{T} \right)}^{-1} A \right)}^{-1} {A}^{T} H b $$


Digital filter design basic principles (IIR/FIR)

Digital filter design is a very large and mature topic and - as you've mentioned in your question - there is a lot of material available. What I want to try here is to get you started and to make the existing material more accessible. Instead of digital filters I should actually be talking about discrete-time filters because I will not consider coefficient and signal quantization here. You know already about FIR and IIR filters, and you also know some filter structures like DF I and II. Nevertheless, let me start with some basics:

A non-recursive linear time-invariant (LTI) filter can be described by the following difference equation

$$y(n)=h_0x(n)+h_1x(n-1)+\ldots +h_{N-1}x(n-N+1)=\sum_{k=0}^{N-1}h_kx(n-k)\tag{1}$$

where $y(n)$ is the output sequence, $x(n)$ is the input sequence, $n$ is the time index, $h_k$ are the filter coefficients, and $N$ is the filter length (the number of taps). The filter taps $h_k$ are also the impulse response of the filter because if the input signal is an impulse, i.e. $x(n)=\delta(n)$, then $y(n)=h_n$ (if the filter's memory has been initialized with zeros). Equation (1) describes a linear time-invariant finite impulse response (FIR) system. The sum on the right-hand side of (1) is a convolution sum, i.e. the output signal is obtained by convolving the input signal with the impulse response. This is always true, but for IIR filters we cannot explicitly compute the convolution sum because the impulse response is infinitely long, i.e. there are infinitely many coefficients $h_k$. One important advantage of FIR filters is that they are always stable, i.e. for a bounded input sequence, the output sequence is always bounded. Another advantage is that FIR filters can always be realized with an exactly linear phase, i.e. they will not add any phase distortion apart from a pure delay. Furthermore, the design problem is usually easier, as we will see later.

A recursive LTI filter is described by the following difference equation:

$$y(n)=b_0x(n)+b_1x(n-1)+\ldots+b_Mx(n-M)-\\
-a_1y(n-1)-\ldots-a_Ny(n-N)\tag{2}$$

Equation (2) shows that the output is not only composed of weighted and delayed input samples, but also of weighted past output samples. In general, the impulse response of such a system is infinitely long, i.e. the corresponding system is an IIR system. However, there are special cases of recursive filters with a finite impulse response. Note that the impulse response is not anymore given by either the coefficients $b_k$ or $a_k$ as in the case of FIR filters. One advantage of IIR filters is that steep filters with high stopband attenuation can be realized with much fewer coefficients (and delays) than in the FIR case, i.e. they are computationally more efficient. However, one needs to be careful with the choice of the coefficients $a_k$ because IIR filter can be unstable, i.e. their output sequence can be unbounded, even with a bounded input sequence.

Filters can be designed according to specifications either in the time (sample) domain or in the frequency domain, or both. Since you've mentioned a cut-off frequency in your question, I assume you're more interested in specifications in the frequency domain. In this case you need to have a look at the frequency responses of FIR and IIR systems. The frequency response of a system is the Fourier transform of its impulse response, assuming that it exists (which is the case for stable systems). The frequency response of an FIR filter is

$$H(e^{j\theta})=\sum_{k=0}^{N-1}h_ke^{-jk\theta}\tag{3}$$

where $\theta$ is the discrete-time frequency variable:

$$\theta=\frac{2\pi f}{f_s}$$

with the actual frequency $f$ and the sampling frequency $f_s$. From (3) you can see that approximating a desired frequency response by an FIR system is basically a problem of polynomial approximation. For recursive systems we have

$$H(e^{j\theta})=\frac{\sum_{k=0}^Mb_ke^{-j\theta}}{1+\sum_{k=1}^Na_ke^{-j\theta}}\tag{4}$$

and you get a rational approximation problem, which is usually much more difficult than the polynomial approximation problem in the case of FIR filters. From (3) and (4) you can see that the frequency response of an FIR filter is of course only a special case of the response of a recursive filter with coefficients $a_k=0$, $k=1,\dots,N$.

Let's now take a quick look at filter design methods. For FIR filters you could take an inverse Fourier transform of the desired frequency response to get the impulse response of the filter, which directly corresponds to the filter coefficients. Since you approximate the desired response by a finite length impulse response you should apply a smooth window to the obtained impulse response to minimize oscillations in the actual frequency response due to Gibbs' phenomenon. This method is called frequency-sampling method.

For simple standard filters like ideal lowpass, highpass, bandpass or bandstop filters (and a few others), you could even analytically calculate the exact impulse response by taking the inverse Fourier transform of the ideal desired response:

$$h_k=\frac{1}{2\pi}\int_{-\pi}^{\pi}H(e^{j\theta})e^{jk\theta}d\theta$$

This integral is easy to evaluate for piecewise constant desired responses, as is the case for ideal frequency-selective filters. This will give you an infinitely long, non-causal impulse response, which needs to be windowed and shifted to make it finite and causal. This method is know as window-design.

There are of course  many other FIR filter design methods. One important numerical method is the famous Parks-McClellan exchange algorithm which designs optimal filters with constant passband and stopband ripples. It is a numerical approximation method and there are many software implementations available, e.g. in Matlab and Octave.

The most common IIR design method for frequency selective filters is the bilinear transformation method. This method simply uses analytical formulas for the design of optimal analog filters (such as Butterworth, Chebyshev, Cauer/elliptic, and Bessel filters), and transforms them to the discrete-time domain by applying a bilinear transformation to the complex variable $s$ (analog domain) which maps the (imaginary) frequency axis of the complex $s$-plane to the unit circle in the complex $z$-plane (discrete-time domain). Don't worry if you do not yet know much about complex transfer functions in the analog or discrete-time domain because there are good implementations available of the bilinear transform method, e.g. in Matlab or Octave.

There are of course many more interesting and useful methods, depending on the type of specifications you have, but I hope that this will get you started and will make any material you come across more understandable. A very good (and free) book covering some basic filter design methods (and a lot more) is Intoduction to Signal Processing by Orfanidis. You can find several design examples there. Another great classic book is Digital Filter Design by Parks and Burrus.


Question about vanishing moments in wavelet transforms

The initial motivation of a sparse representation for a class of signals or images is to pre-separate components with different morphological behaviors.

A traditional model for images is the augmented cartoon: piece-wise smooth + contours + geometrical textures + "noise (unmodelled)".

The mother wavelets (gradient-like operators) can capture isolated singularities at different scale with a few coefficients (and if they oscillate, a bit of textures as well). If the singularities are mixed with a polynomial trend, a standard gradient may also produce coefficients. With a sufficient number of vanishing moments, wavelets somehow "pre-separate" polynomials and singularities, and the polynomial part is cast, subsampled, to a small set of approximmation coefficients.

All in all, you can get few coefficients for singularities, and few coefficients for polynomials between singularities (a long polynomial can be encoded into a few polynomial factors). But still you preserve all the information, especially about the gradients you want to keep.

Wavelets are, under some theoretical conditions, quite optimal for such 1D signals. In practice on 2D data, other aspects come into play and counterbalance the need for vanishing moments: the nature of singularities in images, their finite size, the other wavelet properties, etc.

Resultingly, 1D wavelets work ok for not-too-complex images, but become limited with complicated data, complicated processing, or  high noise conditions.

My experience is that, unless you have specific needs, 2-4 moments suffice in most of the cases, meaning that if finally you use a wavelet with more moments, it will be a by-product of other wavelet properties, and not a strong requirement. There are many options in finding the optimal multiscale  geometric representation.


Simulate butterworth filter with initial condition of state vector in R

If you want to filter signals with an IIR filter, it's best to use the "standard" filter function of the language. That would be filter() for Matlab, lfilter() for Python and for R I believe this is signal::filter().
All of these functions allow you to specify the initial state. The exact details of how this is done, depends on the specific implementation of the filter.
If you don't specify any initial state the functions typically assume it's zero.
If I understand you correctly, you want to initialize the filter with the condition of $s[n] = s[0], n < 0$, i.e. a non-zero constant.
Apparently signal::filter() allows you to specify the arguments. initX and initY. I would assume the length of this arguments must match the filter order. You can set initX= s[0], for a first order filter. For an Nth order filter you would have to repeat s[0] N times. You would also have to calculate the steady state output for a constant input, but since it's a high pass filter, that should just be zero, so initY = 0 shuld work.
CAVEAT: I don't speak R and don't have an R development system to work with. This answer is derived from "best practices" and reading some documentation.


The plot is not of power.  The units of the $y$ axis are dB/Hz. The dB part means that a decibel scale is used, which means it is a scale relative to a reference value.

All that negative values on a decibel scale mean is that the actual power is less than the reference power (because the ratio is less than one).

The higher the decibel value, the more powerful the signal.

Your signal has a major peak at about 4Hz, and lesser peaks at 28Hz and 39Hz.  There might be something happening at about 12Hz, but it's hard to say without more information (or data).


The bandwidth of the RF pulse does not change at all - it is the center-frequency of the RF-pulse that needs to be adjusted. The slice thickness is related to the frequency spread caused by your slice-selection gradient (as you see in the first formula). Hence, to achieve the same slice thickness for each slice, the RF bandwith may not change, but the RF center frequency.


It is common for low-cost direct-conversion (zero-IF) RF receivers to have a DC offset in their A/D hardware, which would correspond to a large spike at the frequency that you're tuned to. That's probably what you're seeing in the waterfall display, as some models of RTL2832U receivers use zero-IF architectures.

For what it's worth, for such a cheap device, I'd say the above is actually pretty good from a DC-offset perspective. You can work around it by tuning slightly away from the signal that you want, then apply the residual frequency offset yourself in your signal processing (it looks like that's what your SDR software package is doing).


You got the shift a bit wrong – a shift in time domain (which you want) is a multiplication with a sinusoid in frequency domain.

This is simply a consequence of the convolution theorem: you want to convolve with a dirac $\delta(t-\Delta t)$ in time, so you need to multiply with the dirac's Fourier Transform in frequency domain.

Thus, what you'd need is

DFT -> $\cdot e^{j2\pi \tilde{\Delta t} f}$ -> IDFT.

Yet, you incur other problems here:

To do a DFT of your whole data, you'd need to first put your whole signal in a buffer. In reality, that means you can't do that to work with live signals.
The DFT will lead to a cyclic temporal shift. I.e. the (fractional) samples that you "shift out" at the end will appear at the start of your output! That means you can't simply divide your signal into smaller DFT chunks and shift each one of them separately.

If you need to avoid that: Use the fast convolution; overlap-add and overlap-save are the two algorithms you should be looking at. It's what I did in a lot of places, when I needed multiple variable delays.

Generally, any (time-)symmetric linear filter has constant group delay – you don't even have to go through frequency domain if you don't like to.


unsigned audio formats, does DC offset matter?

Consider a speaker which wobbles in and out as its driven by a fluctuating audio curve signal ... when the given audio is silent the speaker driver is at rest, not moving ...  that is if the input audio is properly normalized such that it has no DC offset ...  speakers are engineered to deliver optimum fidelity when given zero DC offset audio curve

Now when the speaker is driven by audio which does contain a DC offset the speaker driver is no longer fluctuating in and out centered at its silent resting driver code position ... so when audio is silent yet with a DC offset I would say the speaker driver is held at a physical offset away from its resting position ... at the very least this will reduce the efficiency of the speakers so perceived volume will be lower and I dare say with reduced fidelity

UPDATE googling :  if audio has a DC component offset will silence make speaker cone stay off its resting position ... wisdom is its bad for speakers  https://electronics.stackexchange.com/questions/233230/why-are-dc-signals-bad-for-loud-speakers


What is the name of this digital low pass filter?

I'm new, so I can't add this comment to Matt L.'s answer.

It is not an exponential filter, the equation is actually:

$$ y[n] \ = \  \alpha \, x[n] \ + \ ( 1 - \alpha ) \, x[n-1] $$

So it is a very short FIR filter, not an IIR filter.  I'm not expert enough to know a specific name.

Ced

=============================================
Followup:

I want to thank everyone for their upvotes.  Yes, I've contributed to comp.dsp and I am a blogger at dsprelated.com.

Like all of you, I suspect that the intent of the function was to be an exponential filter and was coded erroneously.

If I were to name the filter as coded, I would call it a "Linear Interpolation Filter", or perhaps a "Subsamplesize Time Shift filter" when $\alpha$ is between zero and one.  It would make more sense as such if the $\alpha$ and $(1-\alpha)$ coefficients were reversed.


Correct way to simulate Rayleigh fading for indoor environments

A discrete-time baseband model of a multi-propagation channel can be written as

$$y[m] = \sum_l h_l[m] x[m-l] + w[m]$$

where $l$ is the index of channel taps if the channel is modeled as a FIR $\{h_l, 0\leq l < L-1\}$. $w[m]$ is AWGN sample.

The channel tap $h_l[m] = \sum_i \alpha_i(m/W) \times \mathrm{sinc}(l-W\times\tau_i(m/W))$ where $i$ is the index of physical propagation path; $\alpha_i$ and $\tau_i$ are the gain and the delay of path $i$, respectively; $W$ is the bandwidth of signal $x(t)$.

The philosophy of this modeling is that due to the decay of $\mathrm{sinc}()$ function, $h_l[m]$ has significant contribution of the path $i$ having delay $\tau_i \in [l/W-1/2W,l/W+1/2W]$. Next the reason of Rayleigh modeling is that these number of that kind of path $i$ for $h_l$ is quite high and by the central limit theorem, $h_l$ is complex Gaussian.

The Rayleigh fading assumption is reasonable for the indoor environment because of the high number of scattering object.

The fact that you have used element-wise multiplication between s_l and rayleigh_h implies that your channel is one tap, i.e. $L=1$. One reason is that bandwidth $W \ll \textrm{delay spread } \tau_m$. Check it.

To model the fact that channel varies in time, you need to consider the notion of coherence time. The channel model you are trying to use is called Rayleigh Block fading. For a quick and dirty approach, just calculate the coherence time $T_c$ (based on the carrier frequency and the velocity), then divide it by the symbol time $1/W$ for the number of symbol that channel is assumed constant. I mean something like

nrep = OSR * Tc * W;
rayleigh_h = 1/sqrt(2)*[randn(1, length(s_l) / nrep) + 1i * randn(1, length(s_l) / nrep)];
rayleigh_h = repelem(rayleigh_h, nrep);
s_l = s_l .* rayleigh_h;

For more sophisticated model, I suggest using the rayleighchan object.


Here's a general example on how to perform the convolution + decimation for an ASK signal.

amplitudes = [1, 0 0 1 0 1 1 0 1 0];
N = 100; % the symbol length.
pulse = 1.2*ones(N, 1); % the Rx pulse (used for correlation)

% the tx pulse (used to recreate your RX signal);
tx_pulse = pulse + 0.2*sin(2*pi*10/N*(1:N)');

% artificially create a signal
signal = zeros(N*length(amplitudes),1);
for i=1:length(amplitudes)
    signal((i-1)*N+(1:N)) = amplitudes(i) * tx_pulse;
end
signal = signal + 0.01*randn(size(signal));

subplot(3,1,1);
plot(signal);
title('RX signal');

subplot(3,1,2);
after_conv = conv(signal, pulse);
plot(after_conv);
title('After convolution');

subplot(3,1,3);
stem(after_conv(1:N:end));
title('After decimation');


It's the correlation coefficient that you get for maximum matching, which is the squared Euclidean 2 norm of the signal $s(t)$, ie. $||s||_2^2$; when you note down the crosscorrelation in terms of integrals, you'll quickly notice that, if you use both signals for correlation, the value for $\tau=0$ is the same as the energy in the signal.


DFT of discrete signals, why do we only analyze frequency bins equal to number of input samples?

A DFT is an invertable (square matrix) transform, thus the complex result vector is the same size as the complex input vector.  But you can analyze as many (additional) frequencies as you want by interpolating the results after the DFT (which requires additional computation that doesn't really increase separation resolution, but can produce a smoother looking plot).

You can analyze less if you want by using the Goertzel algorithm, which can be used to compute any single frequency result bin of a DFT.  (But calculating less than N orthogonal bins won't be invertable or "complete".  And any more than log(N) bins will likely be a lot slower than just doing a full FFT and throwing away what you don't need.)


FFT algorithm gives wrong signs on my output numbers

Reversing the order of the input sequence would provide the same result as what the OP achieved, but not to say this is exactly why it is occurring here. Below shows the details of the DIT algorithm after each stage; comparing each element step by step with the code should reveal the actual error.

Interpreting the OP's results (I assume [x] represents the index) reveals the error is in the core 2 pt DFT butterfly operation.
While the 2 pt DFT is given as
$Y[0] = X[0] + X[1]$
$Y[1] = X[0] - X[1]$
While the OP is doing the following:
$Y[0] = X[0] + X[1]$
$Y[1] = X[1] - X[0]$
With $X[0]$ and $X[1]$ as the 2-Pt DFT input and $Y[0]$ and $Y[1]$ As the 2-Pt DFT output.
  First stage:
  [0] + (1+0i)[2]   = 0+2
  [2] +(-1+0i)[0]   = 2-0   (should be 0-2= -2!)
  [1] + (1+0i)[3]   = 1+3
  [3] +(-1+0i)[1]   = 3-1  (should be 1 -3= -2!)
  Second stage:
  [0] + (1+0i)[2]
  [1] +  (0-i)[3]
  [2] +(-1+0i)[0]
  [3] +  (0+i)[1]


Correct transfer function with 1st order IIR filter

If the formula becomes the same, the response becomes the same – there's no difference to be expected.


Since you don't say what response you're actually getting from Matlab, this is all guesswork.

So let's do the basics: A little formula juggling.

$$\begin{align}A_0(z)&=1\\
A_1(z)&=\frac{-a+z^{-1}}{1-az^{-1}}\\[4em]
H_{LP}(z)&=A_0(z)+A_1(z)\\
   &= 1 + \frac{-a+z^{-1}}{1-az^{-1}}\\
   &= \frac{{1-az^{-1}}-a+z^{-1}}{1-az^{-1}}\\
   &= \frac{(1-a)\quad +\quad (1-a)z^{-1}}{1-az^{-1}}\\[4em]
H_{HP}(z)&=A_0(z)-A_1(z)\\
   &= 1 - \frac{-a+z^{-1}}{1-az^{-1}}\\
   &= \frac{{1-az^{-1}}+a-z^{-1}}{1-az^{-1}}\\
   &= \frac{(1+a)\quad +\quad (-1-a)z^{-1}}{1-az^{-1}}
\end{align}$$

So, taken from these formulas, the filter coefficients for $H_{LP}$ are $b=[1-a,1-a]$, $a=[1,-a]$, and  for $H_{HP}$ it's $b=[1+a,-1-a]$, $a=[1,-a]$ for the canonical form of recursive filters

$$H= \frac{\sum\limits_{i=0}^N b_i z^{-i}}{\sum\limits_{i=0}^M a_i z^{-i}}\text.$$

Plugging that into scipy's scipy.signal.freqz for different $a$ yielded (abscissa: normalized frequency $\in[0,\pi[$, ordinate: magnitude):

$H_{LP,a}$ for $a\in\{0,0.1,\dots,0.9,1\}$:

$H_{HP,a}$ for $ a\in\{0,0.1,\dots,0.9,1\}$:

which both very nicely match the terms low and high pass filters, respectively.


Why is the frequency measured in radians in filter design?

The value of $2\pi$ radians corresponds to the sampling frequency. So the normalized frequency in radians is given by

$$\omega=\frac{2\pi f}{f_s}\tag{1}$$

where $f$ is the actual frequency in Hertz, and $f_s$ is the sampling frequency in Hertz.


Additional prefix and suffix required when windowing an OFDM symbol?

Yes, it is necessary to add an additional pre/suffix to the OFDM symbol that corresponds to the window length when windowing is applied. As the length of the already existing cyclic prefix is usually chosen as the maximum delay spread of the channel no further Inter-symbol interference (ISI) must be introduced or otherwise orthogonality is lost. The windowing method described in the link you referenced is introducing additional ISI by overlapping consecutive symbols.

Furthermore, weighting some samples of the OFDM time-domain signal with different factors is a distortion that also distorts the constellation in frequency domain, as you say.

Note that windowing is not filtering. Windowing is the multiplication of the OFDM signal $x_k$ with some window function $w_k$. In contrast, a filter has a memory.


Wavelet Fidelity Filter (video compression)

Eventually, I found the answer in this white paper. I quote:

  The Fidelity filter was specifically developed for VC-2 using an
  exhaustive search process to select a filter set with simple integer
  filter coefficients, limited to 8-bit resolution, and offering good
  alias rejection in both analysis and synthesis operations.


Oversample or undersample correlated noise?

I would strongly suggest oversampling at a frequency much higher than the bandwidth of concern, and specifically setting up the test with proper filtering to ensure that there is insignificant energy beyond half the sampling rate or $f_s/2$ (assuming real sampling). This would imply using a low pass filter prior to the sampler and also reviewing the resulting power spectral density that the energy is sufficiently lower as you approach the Nyquist frequency of $f_s/2$.  The analog input frequency that extends to $f_s/2$ is referred to as the first Nyquist band.

You may be "bandpass" sampling in which case "undersampling" would apply with the same guidance above to "oversample" at more than 2x the bandwidth of interest with sufficient margin for filting. In this case the filter is a bandpass filter centered on the higher Nyquist zone (the first Nyquist zone for real signals is $f=0$ to $f_s/2$ as noted above, the second Nyquist zone is $f_s/2$ to $f_s$, the third Nyquist zone is $f_s$ to $3f_s/2$ etc all of which may be undersampled.

If this is still confusing, please see my post here with further details:

How do you simultaneously undersample and oversample?

However in all cases it is imperative to be filtering the signal of interest and to oversample it with a sampling rate > 2x the resulting bandwidth of the filter.
The reason is when you undersample below that bandwidth, while you will still be able to measure the overall total integrated power density, you will have lost all information as to what frequency your noise issues are at which would typically be of great interest in dealing with the noise. When sampling, everything that is in a band over a multiple of the sampling rate will fold into the same frequency band in your digital spectrum that extends from $-f_s/2$ to $+f_s/2$ (which is equivalently $f=0$ to $f_s/2$ for real signals).  For example, every frequency $Nf_s + \Delta$ will alias to $\Delta$ in your digital spectrum.

The easy way to see that the total integrated power density is conserved is to consider sampling white Gaussian noise and take the histogram of the results to visually see the bell curve and related standard deviation, the square of which is the variance proportional to the total power in the signal. From this visual it is easy to convince yourself that you can sample at ANY rate and with enough samples the histogram would have the same shape, standard deviation and variance. (This is true as long as the sampling clock is uncorrelated to the signal which would be the case for white noise.) But the content of the signal itself is all the energy at the input to the sampler at all frequencies where energy exists, so the actual analog frequency is not recoverable if not all isolated to one Nyquist band.


Sketching Phase Spectra Using Group Delay and Magnitude Spectra Informations

If the magnitude spectrum is symmetric

$$M(\omega)=M(-\omega)\tag{1}$$

(as I assume), then your system is real-valued. The phase response of a real-valued system is asymmetric:

$$\phi(\omega)=-\phi(-\omega)\quad(\mod 2\pi)\tag{2}$$

This means that there can be two cases:

The phase goes through zero at $\omega=0$, i.e. the phase is given by $\phi(\omega)=-c\cdot\omega$, where $c$ is the constant group delay.
The phase is $\pm \pi$ at $\omega=0$, which means $\phi(\omega)=-c\cdot\omega\pm\pi$.
The phase jumps at $\omega=0$. This is only possible if the magnitude has a zero at $\omega=0$, i.e. $M(0)=0$. The phase jumps by $\pi$, which is simply a sign change of the (bipolar) amplitude function. In this case the phase is given by $$\phi(\omega)=\begin{cases}-\pi/2-c\cdot\omega,&\quad\omega>0\\\pi/2-c\cdot\omega,&\quad\omega<0\end{cases}$$


(Rough) camera calibration of manipulated images

1) It seems to me that you don't care about solving for correspondences or obtaining any high level paradigm in that sense. So I would treat this problem as a minimal problem and redirect you to check:

http://cmp.felk.cvut.cz/minimal/

Maybe specifically, 4-point method with unknown focal length:
http://cmp.felk.cvut.cz/minimal/p4pfr.php

You could find many variants in that page.

2) Well, this problem involves many issues of multiple view geometry theory. However, in this case, an uncalibrated approach would be more suitable.

You could obtain point correspondences from SIFT-like features. As your images are downloaded from social platforms, it is likely that they will contain a lot of such feature points. Then, you could match the descriptors in a RANSAC-like fashion and retain the correct matches. From there it is possible to compute $F$ the fundamental matrix. With further assumptions, constraints or knowledge about the scene you could compute the scale as well as  the 3D point coordintates. Some examples here:

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.31.4741&rep=rep1&type=pdf

Under certain constraints it is also possible to recover the camera matrices:

http://users.cecs.anu.edu.au/~hartley/Papers/eccv92/Higgins/higgins.pdf

After all, I would recommend you to read:

http://www.robots.ox.ac.uk/~vgg/hzbook/


Speech recognition and comparison algorithms

What you want to do sounds like standard speech recognition. I recommend you check out the Open Source JAVA library from CMU, Sphinx4. Run the Transcriber or Dialog demos and you can see it in action doing some basic recognition of digits & simple english commands from a small grammar. Note that many of your algorithms listed above, fit into different parts of a speech recognition system (frontend processing, search manager, acoustic scorer, etc)

Check out this paper for a good overview (and to see how all your algorithms above relate)


What mother wavelet should be chosen as an alternative to the STFT?

Just because the math is easier, I might recommend the Morlet or possibly the Chirplet as the mother wavelet.


Here's an attempt to answer your first question:

  "An example of a 94% overlap transform processing could be demonstrated by a 279 sample frame with 262 samples overlapped from each frame into the next.

A spectrogram is really just a sequence of spectra.  In this case, each individual spectrum is taken of 279 time samples (which I will call a frame).  Each subsequent frame consists of the last 262 samples of the previous frame, plus the next 17 time samples.

That means percentage of common samples between subsequent frames is:

$$
\frac{262}{279} = 93.9068 \%
$$

  At 12.5 MHz sample rate, each spectrum frame would occupy 23.76 µs.

We have 279 samples in each frame.  The distance between time samples is

$$
\Delta t = \frac{1}{12.5 \times 10^6} = 80 {\rm ns}
$$

Add 279 of those together and we get

$$
279 \times 80 {\rm ns} = 22.32 {\rm \mu s}
$$

It looks like someone can't multiply:  the number 23.76 µs is actually $297 \times 80 {\rm ns}$.

  Each spectrum would start 1.36µs after the previous one...

Each frame starts 17 samples after the last.  This is $17 \times 80 {\rm ns} = 1.36 {\rm \mu s}$

  ... and could show a 10 MHz span." Could anybody please help to explain what is 10 MHz span and why why each spectrum could show a 10 MHz span?

It's not clear what they are meaning by this.  If you sample at 12.5MHz, then the total possible frequency span is 12.5MHz (if you include positive and negative frequencies).  There may be some allowance for guard bands that reduces the usable bandwidth from the ideal.


Does the sound volume affects the MP3 file size?

Depends on the details, but typically not. For the sake of discussion, we assume CD quality for lossless, i.e. 16 bit, 44.1 kHz.

Reducing volume may or may not reduce information. It depends on what number format it's done in. If you convert 16-bit input to, to for example, 32 floating point, no information would be lost. It may still be lost later, depending on what happens to the data down the road
If you do it directly to the 16-bit input data, you have to
basically throw away the LSB (least significant bit), and your
bit depth has been reduced to 15 bit.
This has little or no effect on MP3 encoding. Most encoders start with a target bit rate and then minimize the encoding artifacts within these constraints. If you encode 10s wave file at 320kb/s, the resulting file will be very 3.2 Mbit long. The thing that varies with content is the sound quality and not the file length.
The type of "information" that encoders use are "masking" and "signal to noise ratio". These are largely independent from the overall scale, so shifting volume doesn't affect the encoders much.
Most encoders these days come in two variants: constant bit rate and variable bit rate. With "variable bit rate" you may see the file length fluctuate a bit, but generally not by much. With constant bit rate, the file lengths will always be constant.


Thoretically, what makes multi-pitch detection so difficult?

First of all:

  mostly the frequency with the largest magnitude is desired frequency

This is a very naive and flawed approach, as what you are measuring here is very different from the perceptual notion of pitch. Not to mention that you hit the limited resolution of FFT...

Some of the challenges one has to face in multi-pitch detection:

Octaves are difficult to tell apart. Play two Cs one octave apart, and the harmonics of the highest note reinforces the harmonics of the lowest note, rather than create new peaks in the spectrum, to the point that one could hear one single "C" note with a different timbre. When mixing several notes together, some of their partials can reinforce or cancel each other.
Pitch detection can no longer be seen as a problem of period estimation. The most successful monophonic pitch detection algorithms are time-domain and based on metrics like the average magnitude (or square) difference function, and the underlying assumption is that the signal has a salient periodicity - it is almost similar to itself when translated in time. These methods cannot be applied to polyphonic signals, as the period of a mixture of periodic signals is the LCD of the components' periods. Knowing this LCD does not give away the factors...
The fact that the number of notes is not known in advance increases the size of the search space, and makes it an ill-posed problem. You are no longer ranking hypothesis like "this is a C3" vs "this is a C#3", but instead, "this is a C3 and a G4" vs "this is just a C3" vs "this is a C2 and a C3". For example, take a single piano note. One could describe its spectra as a "single C3 piano note". Or one could describe its spectra as "a C3, a C4, a G4, etc... all played together by an instrument producing purely sinusoidal notes". Distinguishing among these two options requires prior knowledge of the spectral shape of a piano note.


Sound behaves (more or less) like a wave. When the positive peak hits the membrane it is moved in one direction and afterwards "sucked" in the other. Just as waves on the beach do.
Each positive peak generates a equal negative peak (in a perfect world...). Superimposition of different waves makes it hard to see this directly.

The average is not necessarily the atmospheric pressure, but is a combination of it in combination with the property of the electrical circuit of the microphone.

For your fft you have to consider positive as well as negative values. However, you may want to get rid of the average value if you do not care for very low frequencies.


M-point FFT Amplitude of N-samples of a sinusoid

In general the FFT of a sinusoid function is NOT a delta function. The FFT's maximum will be around the sine frequency but the exact value will be dependent on how well the frequency of the sine lines up with the FFT grid and there will be non-zero values for most other frequencies.

You only get a delta function if the sine frequency is an integer multiple of your FFT bin spacing (= sample rate divided by FFT length). In this case the amplitude will N/2 (for conventional FFT scaling) for both the positive and the negative bin.

If your FFT length is shorter than the sine wave you have to truncate and the above considerations apply. If the FFT length is longer, you will probably zero pad. This will often result in substantial increase of the energy at the "other" frequencies and the effects depend heavily on the details of the zero padding.

In all cases Parseval's Theorem holds: total energy in the time domain equals total energy in the frequency domain. That's how you can derive an amplitude of N/2 for the special case that the sine wave frequency falls on an FFT bin: A sine wave of amplitude 1 has an RMS value of sqrt(.5). A buffer of N samples has a total energy of sqrt(.5)*N. This is the total energy. In the frequency domain this gets split between the bins at +f and at -f. Each bin receives have of the total energy which is sqrt(.5)*sqrt(.5)*N = N/2.


Compare Overlapping FFTs to a longer FFT

That depends on your application, signal and noise.
Increasing the FFT length means better frequency resolution.
It can be shown that the FFT of a zero-mean white noise is also white noise.
Now, given that observation, this is the motivation to performing multiple shorter FFTs since the noise averages and its variance decreases by a factor of the number of FFTs.
For a stationary signal, e.g. a sine wave, the FFT should not (theoretically) change over time. Thus, averaging would have no effect. So in this case, since the noise decreases, you should expect a gain in SNR by a factor of the number of FFTs.

However, one should also take into account the scalloping loss (or spectral leakage) which can reduce the SNR significantly, and choose the FFT length correctly, according to the application.

Personally I always use overlapping if possible.


MLE formulation -- confusion regarding the terms in the equation (Part1)

Hi: I don't have time to carefully look at everything right now and I'm unfamiliar with circular complex random variables. But, if it's similar to regular normal rv's and the two variables are independent , then their sum has mean zero and their variance is $ 2 \sigma^2 + 2 \sigma^2_{v}$. So, the likelihood  can be written having mean zero and that variance.  But it's not clear to me what you're estimating. There aren't any coefficients in the model ( implied one's I guess somewhere or maybe you are doing bayesian ? ) so what is actually observed ? $x_{n}$ or
 $v_{n}$. That has to be known if you want to write the likelihood correctly.
If $x$ is viewed as fixed and known and not a random variable, then what you wrote is missing a factor of 2 in the factor before the exp. The $2 \pi$ has to be multiplied by the variance which is $2 \sigma^2_{v}$. Also, it's not clear to me how X can be viewed as fixed since you said in the question that it was an RV also. If it is an RV and also not observed, then that has to be accounted for in the likelihood.  That's all I can tell you and even what I'm saying may not be correct if complex means something totally different. Good luck.


Doubling the frequency of a sound fragment in real time?

After loads of research, I found a really good article about Short Time Fourier Transform (STFT), which can apply FT to a really small buffer, change it's frequency and apply IFT. There's even an implementation in C. You can find it here : http://blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/

As a "side effect", you also know the pitch of your signal in almost real time

HTH


once again,confusion between phase and group delay

If you are looking for a frequency-independent delay applied to any given input signal by the filter (apart from amplifying and attenuating certain frequency components), then you won't be able to find it because there is no such delay. As you can see in your plots, group delay and phase delay are generally frequency dependent. Furthermore, for general input signals, the terms group delay and phase delay have no meaning.

There are two cases for which the terms group delay and phase delay have a clear meaning. The first is for narrow-band input signals

$$x(t)=a(t) \cos(\omega_0 t)$$

where $a(t)$ is a  low-pass signal. If the system has an approximately constant amplitude response in the frequency range of $x(t)$ (i.e., around $\omega_0$) which is $\left|H(j \omega_0)\right|$, and if its phase is approximately linear in this frequency range, then it can be shown that the output signal is approximately given by

$$y(t) \approx \left| H(j \omega_0) \right| a(t-\tau_g(\omega_0)) \cos\left( \omega_0(t-\tau_\phi(\omega_0)) \right)\tag{1}$$

where $\tau_g(\omega_0)$ is the system's group delay evaluated at the carrier frequency $\omega_0$, and $\tau_\phi(\omega_0)$ is the phase delay at $\omega_0$. So in this case the group delay is the delay of the envelope, whereas the phase delay equals the delay of the carrier. A special example of such a narrow band input signal is a pure sinusoid. From (1) it is clear that the delay of a pure sinusoid when passing through a linear time-invariant (LTI) filter equals the phase delay evaluated at the sinusoid's frequency (and not the group delay). Note that in the case of a pure sinusoid (i.e. $a(t)=\mbox{const}$), Equation (1) becomes exact for any LTI system.

The other case where these terms are meaningful is when the system has exactly linear phase. This is only possible for FIR systems. For FIR filters with even symmetry of the impulse response, the group delay and the phase delay are equal and constant. In this case, they simply specify the delay of the input signal. This is possible because both group delay and phase delay are constant, i.e. frequency independent.

From your question, I believe that you are looking for such a delay. However, as explained above, such delay only exists for linear-phase FIR systems, not for the Butterworth filter you are considering.


What is the best input for de-noising autoencoder for sound data?

As you mentioned MFCC features are one of the best features to represent audio as it captures both the time and frequency variations in the audio clip.You can get more details about MFCCS features in the below link:
http://practicalcryptography.com/miscellaneous/machine-learning/guide-mel-frequency-cepstral-coefficients-mfccs/

You can import python_speech_features in python and use them to directly find the MFCCs of audio which will be in a matrix form.

You can feed the noisy audio features as input and the clean audio features as output to Deep Denoising AutoEncoder(DDAE) and then you can test on test datasets.

If required you can reconstruct the predicted mfcc output of the DDAE and check the extent to which your architecture works.

You can use Griffin Algorithm for this. Refer the link below for details about the reconstruction.
https://timsainb.github.io/spectrograms-mfccs-and-inversion-in-python.html


Assuming $x(t)$ swings between -1 and 1, $\sin \frac{\pi}{2} x(t)$ might do the trick.


If you are trying to find a signal inside another signal, Correlation is a useful technique.  You will get a large peak anywhere the signals match.  The easiest way to do it is to flip the shorter signal, such that samples 0 & N-1 are swapped, 1 & N-2, etc, then convolve the two signals.  This should attenuate any effects of noise or edge condition resulting from resampling or fades, unless the phase data got mangled, in which case you’ll probably need something more sophisticated.


A 23 kHz sine wave sampled at 48 kHz, then FFT

Since you used 2.3k twice, I am assuming you meant that with the original 23k value.

Converting your frequency:

$$ 2.3k \frac{cycles}{second} \cdot \frac{ 64 \frac{samples}{frame} }{ 48k \frac{samples}{second}  }  \approx 3.07 \frac{cycles}{frame}$$

Because this is not an integer value, the peak DFT value will be at bin 3 and there will be values in all the other bins, which is commonly called "leakage".  Being a real valued signal, the upper half of the DFT will be a mirror image with complex congjugate values, so there will also be a peak at bin 61 (64-3).

Your 100 sample figure is irrelevant if you are only using 64 of them for the DFT.

Signals have to be a whole number of cycles per frame in order for the DFT not to have leakage.  The other answers presume you understand that the DFT is a dirac delta, windowed, convolved version of the continuous FT.  This is not an easy way to understand it, but it is how it is usually taught.  It also leads to the misconception that the leakage is described by what is known as the sinc function $ \frac{\sin(d)}{d} $.  That is not true.  The fall off for a complex valued signal is described by what is known as the aliased sinc function $ \frac{\sin(Nd)}{N\sin(d)} $.  Another name is "Dirichlet function".  A real valued signal is actually the sum of two complex signals, so that function doesn't apply either.

You can find the exact formula for the bin values, the peak and "leakage", in my blog article "DFT Bin Value Formulas for Pure Real Tones".

You can also find a much more intuitive explanation of why leakage occurs in another one of my blog articles "DFT Graphical Interpretation: Centroids of Weighted Roots of Unity".

Hope this helps,

Ced


Help with filter stability, ROC and output of a system

For any filter with AR(1) component the coefficient must be smaller than in order to be stable. Otherwise any input is infinitely multiplied by a non decreasing factor.


How to scale the amplitude of an FFT of a chirp signal

First, thanks for all your inputs. I did some research and I found my mistake. The spectrum of a linear chirp can, according to Wikipedia (search for: "chirp spectrum wiki") be analytical described in the frequency domain. The formula for the amplitude

gave me the hint. I just used the constant factor and a factor of samplingfrequency/numberofsamples to match the FFT. I do not know exactly why, but it seems to fit.

Following Picture shows the result:

Following code can be used to simulate it:

    %% Test script
close all
clear all
clc
% create timevector
L =1000 ;
aux = 0:1:(L)-1;  % indexvector
% get fs
fs = 1000;
% create timevector
t = (1/fs)*aux;
f0 = 1;
f1 = 100;
% chirp in timedomain
SignalChirp = chirp(t,f0,max(t),f1);
% Window length
nfft= length(SignalChirp);
%generate the vector of frequencies
halfn = floor(nfft / 2)+1;
deltaf = 1 / ( nfft / fs);
ffft = (0:(halfn-1)) * deltaf;
% perform FFT
X = fft(SignalChirp,nfft);
magni1(1) = abs(X(1)) ./ (nfft);
magni1(2:(halfn-1)) = abs(X(2:(halfn-1))) ./ (nfft / 2);
magni1(halfn) = abs(X(halfn)) ./ (nfft);

figure(1)
subplot(2,1,1)
plot(t,SignalChirp,'k');
title(['Chirp Signal']);
xlabel('Time(s)');
ylabel('Amplitude');
subplot(2,1,2)
plot(ffft,magni1),hold on
duration = max(t);
c = (fs/L);
% https://en.wikipedia.org/wiki/Chirp_spectrum
scalingLinChirp = sqrt(duration/((f1-f0)));
plot(ffft,ones([length(ffft) 1])*scalingLinChirp*c), hold on
xlabel('Freq(Hz)');
legend('FFT Result' ,'Calculated Amplitude' )
ylabel('Amplitude');

If I use this scale factor to my FFT of my accelerometer data it matches the calculated theory much better than before. Sadly I can’t show you guys another picture due my low rep.

Thanks for all your Inputs!


Remove a signal that is known to repeat every 4 samples at a known phase i.e. exactly 500 Hz signal in 2000 Hz data

Currently we would employ a Butterworth bandstop of 490-510Hz but this doesn't completely cut-out 500Hz under some conditions

Use a notch filter instead. It will completely remove the 500Hz component and you can adjust the Q to dial in the width vs latency or group delay distortion

  Also if we make the filter too narrow it has a very long latency.

That's unavoidable: the narrower/steeper a filter is in one domain, the wider it is on the other.

  if a similar/same method is exists to subtract/remove the 500Hz signal from the signal?

Yes, many. For example you can use phase locked loop (PLL) to determine both phase and amplitude of the 500 Hz component and then simply subtract it out of the signal. This will cause NO delay or latency to the remaining signal. If your 500Hz has any type of amplitude or phase drift, the PLL will simply track it and you can adjust the speed of PLL controller to match the drift.

If your 500 Hz noise signal is guaranteed to be EXACTLY a 1/4 of the sample rate, that gets really simple. Just heterodyne it down with a cosine and a sine at 500 Hz, that gives you the real part and the imaginary part. Calculate amplitude and phase, create sine wave and subtract.


Let's consider the complex signal
$$\phi_c(t)=Ae^{j(\omega_ct+km(t))}\tag{1}$$
The power of $(1)$ is readily computed as
$$P_{\phi_c}=\lim_{T\to\infty}\frac{1}{T}\int_{-T/2}^{T/2}\big|\phi_c(t)\big|^2dt=A^2\tag{2}$$
IF $(1)$ were the analytic signal corresponding to the real-valued signal $\phi(t)$ then the power of $\phi(t)$ would simply be $P_{\phi_c}/2=A^2/2$.
In general $(1)$ is not analytic. However, in all practical cases of angle modulation, a sufficiently high carrier frequency makes $(1)$ analytic for all practical purposes, i.e., there are no measurable negative frequency components. Or, equivalently, a practical angle modulated signal is a bandpass signal and has no components at or close to DC. If this is the case, the power of $\phi(t)$ is indeed given by $A^2/2$.
In all other (artificial) cases, the power of $\phi(t)$ cannot generally be shown to equal $A^2/2$.


Butterworth filter for signal processing

First of all when you work with MATLAB I suggest you to use a time axis t.

[x,Fs] = audioread('go10.wav'); % to read the data
t=(0:length(x)-1)/Fs;
plot(t,x) ; % Now when you plot you can see the time in seconds

The code seems to be correct. A butterworth bandpass filter, like others bandpass filter keep only the frequencies in the interval from beginFreq and endFreq. It is impossible to make a perfect filter, so there are a lot of filters like bessel, cebicev, butterworth,... each one have some particular properties and the choice of the better filter depends on the application. In sound processing butterworth filter is quite good. The human voice goes from 100Hz to 1100Hz, if you want to keep the voice in your signal you should change your values to:

beginFreq = 100/ (Fs/2);
endFreq = 1100/ (Fs/2);

A passband filter, like this butterworth, keep only the frequencies in the range you have chosen, so if you have noises in your signal that can't are human voices this filter will remove them. However this filter keeps all noises that have frequency in the vocal range.


Fields of math needed for digital filter design

To get started:

Complex numbers

The frequency response of a filter is easier to understand complex-valued, describing both the magnitude frequency response and the phase frequency response. You will be able to understand poles and zeros, which can be complex. Complex numbers enable you to have negative frequencies, which will make math simpler.

Trigonometry

$\sin$, $\cos$ and their relation to the complex exponential $e^{i\alpha} = \cos(\alpha) + i\sin(\alpha)$ are important. Sinusoidal functions will be passed through filters with only their amplitude and phase affected.

Differentiation

To find at what frequency a simple filter peaks or dips, you can solve at what frequency the derivative of its magnitude frequency response is zero.

Integration

Integration is needed for the Fourier transform and the inverse Fourier transform.

Fourier transform

Fourier transform enables you to go from an impulse response to a frequency response and back. Also things you do in the time domain often have a simple counterpart in the frequency domain, and vice versa.