My 3x3 Weighted Median Filter doesn't seem to be improving my image quality In order to reduce noise in my training dataset,I attempted a WMF whose weights are shown in a 2-D array as follows(values finally get normalized by division by 15:
[1 , 2 , 1
 2 , 3 , 2
 1 , 2 , 1]

My input image is an RGB image size (128,128,3). The pixel values are integers in the range of 0-255.
Isolation of each of my three channels is done as:
X_R=X_dummy[:,:,0]
X_G=X_dummy[:,:,1]
X_B=X_dummy[:,:,2]
R_padded=np.pad(X_R,[mask_start,mask_start],'symmetric')
G_padded=np.pad(X_B,[mask_start,mask_start],'symmetric')
B_padded=np.pad(X_G,[mask_start,mask_start],'symmetric')

Where mask_start=2
Convolution for a single channel is shown as:
for x in range(mask_cent,row_pad-mask_start):#The upper limit thing is not included. So Highest value is row_pad-mask_start-1
    for y in range(mask_cent,col_pad-mask_start):#These loops seem to be fine
      patch_selected_R = R_padded[x - mask_start-1: x + mask_start, y - mask_start-1: y + mask_start]
      a1_R= W*patch_selected_R
      med=np.median(a1_R)
      med=med*15
      b[x,y,0]=med

My images, which are shown below look like a colour correction rather than a denoising:

I am unsure of what is going wrong right now.
Edit: Something went wrong, this is my imgur link
https://imgur.com/gallery/1d4V5vr
Edit 2: new code, its resulting output is available through this link
https://imgur.com/Dr9K1bU
:
a1_R=patch_selected_R
      for i in range(1,2):
        a1_R=np.append(a1_R,patch_selected_R[0,1])
        a1_R=np.append(a1_R,patch_selected_R[1,0])
        a1_R=np.append(a1_R,patch_selected_R[1,2])
        a1_R=np.append(a1_R,patch_selected_R[2,2])
      for i in range(1,3):
        a1_R=np.append(a1_R,patch_selected_R[1,1])
      b[x,y,0]=np.median(a1_R)

Apparently, you intend to weight the patch pixel per pixel in:
a1_R= W*patch_selected_R

and then apply a median on the piece-wise product, and divide by 15. By the way, beware of the * product and the input type of the array, you don't want a matrix product.
This is not the weighted median I know of. To me, the elements in  patch_selected_R should be duplicated, triplicated, $n$-plicated with respect to the corresponding integer weight in the mask. On a smallest example: if $w=[1,3,1]$ and $m=[1 4 3]$, the intermediate (longer, of length $1+3+1=5$) buffer is (before sorting):
$$b=[1,4,4,4,3]$$
thus $[1,3,4,4,4]$ after sorting, and the weighted median is $4$ (instead of $3$ with the classical median).
The resulting array here is indeed 15 times larger, and you then compute the median, and replace the central value, with no need to further divide by 15. I am thinking of an interpretation of your result. Meanwhile, the theoretical origin of the weighted median is describe here:
What is the advantage of weighted median filter over median filter?
Note: there are combinations of linear and nonlinear filters, like the mean-medians filters.


Is it possible to estimate variance of noise for a step answer signal? I know there is not possible to find the true noise of a measured signal. The only way to "find" the noise is to estimate the noise. Noise has the mean 0, but the variance varies.

So assume that we have a signal which look like this:

The clean signal is $y(t)$ and the noise is $e(t)$. In this case, $e(t)$ has the variance 1 because it's generated by Octave/MATLAB function 0.05*randn().

So my question is:
How can i estimate variance of this noise from this signal?

Hi: In order to estimate the variance, you need to have an underlying model for your signal. So, suppose that the model is

$y_{t+1} = y_t + \epsilon_t$ $~\forall ~ t = 1,\ldots n $.

assuming that $E(\epsilon_t) = 0$ and $var(\epsilon_t) = \sigma^2$.

In this case, you would difference your data in order to get estimates
of  $\epsilon_{t}$ at each time $t$, and call the estimates
$\hat{\epsilon}_t, ~\forall i = 2, \ldots n$.

Then, then an unbiased estimate of $\sigma^2 , \hat{\sigma}^2$, $= \frac{1}{n-2} \sum_{i=2}^{n}  (\hat{\epsilon_i} - 0)^2$.

But that estimate is critically dependent on the assumed model for the response $y_{t}$ ( this model is sometimes used in finance for the log price at time t+1 and is referred to as a random walk model ).  A different model will lead to a different estimate of the noise because the assumptions about the noise will be quite different in each case. I hope this helps some.


Implementing a Matlab filter to cancel cross talk between two sensors I have two sensors that have a cross talk between them. I would like to cancel the cross talk. For this I recorded two tracks, where sensor no.1 (called x1) has some input, and sensor no.2 (called x2) is zero.

I was told to make the following operations:

  Since your cross talk is small:

  1) Simply measure the transfer function $H_{ba}(\omega)$ directly.

  2) Subtract it out as follows:

  2.1) Measure the transfer function from signal A to sensor B when signal B is 0.

  2.2) Create a filter from this transfer function (FIR or IIR, depending on it's shape).

  2.3) Now you can measure and subtract a filtered version from sensor signal A from sensor signal B:

  $y_b'(t) = y_b(t)-h_{ab}(t)*y_a(t)$. Where $h_{ab}(t)$ is the impulse
  response of your cross talk filter and $*$ the convolution operator.

Here is a Matlab code that I wrote for this purpose, but its performance is really bad. If someone has a suggestion what I did wrong, or can suggest a new code I will be grateful. If I may say, a friend helped me to implement this, and I think that we did not implement the filter part, but maybe I am wrong.

% The files designated to measure the transfer function.
% x1 has some input. x2 is is zero, and affected by the x1.

x1 = wavread(file1);
x2 = wavread(file2);

%The transfer function.

Hab = fft(x2)./fft(x1);

hab=ifft(Hab);

% Some new tracks, where x3 is recorded from sensor no.1 and x4 from sensor no.2.
% All the files have the same length.

x3 = wavread(file3);
x4 = wavread(file4);

c=conv(hab,x3);

% Subtracted signal

x4 = x4 - c(1:length(x4));

The following code did not subtract the cross talk on x4.
Here are some images explaining the problem. For some frequencies the above code actually did a good job and subtracted the cross talk, and for some it did the opposite and increased the energy. plot(abs(fft(x4))):

Good - Frequency energy subtracted:

Bad - Frequency energy increased (!):

Thank you!

Try using other Blind source  separation(BSS) algorithms,even Kernel BSS could be applied to remove if it is nonlinearly mixed.


Technical Term for the Sum I + Q Most likely a very simple question, but I haven't been able to find a good answer.

When working with analytical signals in general, and software defined radio in particular, a common operation is to take the sum real(z) + imag(z). Is there a technical term for this sum?

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

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

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

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

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

These expression can be derived by applying the trigonometric identities

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

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


Mulitplication in Frequency Domain as a Similarity Measure of Audio Clips I have an audio clip and I copy a 256-point part satarting at 10th second. Then I create a 256-point frame on the original clip. At every itreation, I calculate the DFT of the frame, multiply it with the DFT of the interval I copied, take the mean of the result and append it to a list and shift the frame 1 point right. When I graph the means, I know I should get the maximum value of convolution when the copied interval's DFT is multiplied with itself. And this is supposed to show that convolution can be used as a similarity metric. However, when I apply the explained procedure, I can't get a peak that is distinguishable from any other value. What I have done is as follows:
    import numpy as np
    import scipy.io.wavfile
    import matplotlib.pyplot as plt

    rate1, data1 = scipy.io.wavfile.read('Africa.wav')

    data1 = np.array(data1, dtype=np.float64)

    interval_1 = data1[rate1 * 10: rate1 * 10 + 256]

    dft_1 = np.fft.fft(interval_1)

    cv = []

    for i in range(data1.size - 256):
        dft = np.fft.fft(data1[i: i + 256])
        Y = np.multiply(dft, dft_1)
        Y = np.abs(Y)
        Y = np.mean(Y)
        cv.append(Y)

    cv = np.array(cv)

    plt.figure()
    plt.title("Convolution with 256 point sample at 10th second")
    plt.xlabel("Samples")
    plt.ylabel("Amplitude")
    plt.plot(cv)
    plt.show()

I get the following graph:

The fft product alone is not the convolution, but the frequency domain of the convolution. To complete the operation the OP must also take the inverse FFT to get a circular convolution result.
$$CONV = \text{ifft}(\text{fft}(a) \text{fft}(b))$$
However, similarity would be determined using correlation not convolution. To do this, simply complex conjugate one of the FFT results as follows:
$$XCORR = \text{ifft}(\text{fft}(a) \text{fft}^*(b))$$
The above is the cross-correlation function (using circular convolution). The result is the correlation of $a$ and $b$ at repeated circular shifts in time.


How to extract unknown packet from collision between said packet and another, known packet? I need to apply successive interference cancellation (SIC) to an RF signal which is the result of collisions between 2 or more unsynchronised packets. For this question let's just assume that the number of packets is 2.

In SIC, you store the collided packets until one of the packets are repeated and use the 'clean' packet to extract the other packet from the collided data.

My question is, what techniques can be used to extract an unknown packet from a collision between that packet and a known packet?

How do you deal with the packet sources not being synchronised?

Going through the possible tags for this question, it seems that 'deconvolution' is relevant.

In case it is relevant, the context is RFID tag inventory, the channel access mechanism is framed slotted ALOHA (FSA) random access, the carrier is 868 MHz (EU) and subcarrier is 256 kHz.

Thank you.

As JRE said, deconvolution is not what you need.  You need to subtract the known signal from the combined signal to leave the unknown signal by itself.

If the signals are baseband signals (I'm aware that yours are modulated- I will get to that later in the answer) then you need to know when the known signal started and how strong it is.  You can get an estimate of both of those things through cross-correlation.  You cross-correlate the combined signal with a clean reconstruction of the known signal.  You should get a strong peak at the location of the known signal.  From this you can determine when the signal started.  You can also estimate the signal amplitude using the peak's amplitude.

$$
y[n] = \sum\limits_m x[m]r[n+m]
$$

This is the definition of cross-correlation, where x[n] is the combined signal, r[n] is the reconstructed known signal, and y[n] is the cross-correlation.

$$
y[n] = \sum\limits_m (s_1[m] + s_2[m] + noise[n])r[n+m] \approx \sum\limits_m s_2[m]r[m]
$$

In this equation I have broken x[n] into its constituent parts: $s_1[n]$ (the unknown signal, $s_2[n]$ (the known signal), and noise.  If we assume that $s_1[n]$ does not correlate strongly with $s_2[m]$ then we can simply disregard $s_1[n]$ and the noise term with the understanding that they will both likely add some error to the result.

So, the peak's amplitude is approximately equal to $\sum\limits_m s_2[m]r[n+m]$.  Since you know $r[n]$ completely, you can use this value to determine the amplitude of $s_2[n]$.

Once you know the start time and amplitude of $s_2[n]$, you can subtract it out of $x[n]$, leaving only $s_1[n]$, the noise term, and whatever error was in your estimate of $s_2[n]$.

If the packets are modulated the same basic process is followed, except you will also have to find/determine any carrier offset in $s_2[n]$, and it's phase.  The phase you can get from the cross-correlation peak's phase.  The carrier offset you can find using the fourth power technique (for PSK signals), or by searching for the best cross-correlation result at various carrier offsets.


understanding time constant meaning in signal processing from my DSP book,i am reading that  exponential signal $x(t)=e^{-a*t}$ has some time constant value $c$ let say if $e^{-a*t}=e^{-t/c}$ which means that $a=1/c$,but i want to understand   what is general formula to find  time constant for other signals and also what does it express?i meant what is real mining of time constant in signal?clearly it cant be value where signal repeats itself,because  in this case we call this number period,thanks in advance

There is no general formula to calculate time constants for all signals. Time constants are defined for some signals. For your exponential signal

\begin{equation}
e^{-a \cdot t} = e^{-t / c}
\end{equation}

it just means that always after the time $c$ has passed, the value of the function has been reduced by $1 / e$. A more common interpretation is half-life defined as

$$t_{1/2} = \frac{ln(2)}{a}$$

which is the time after the function has gone from 1 (at time zero) to 1/2 (at time $t_{1/2}$) and so on.

However, there is no general time constant that can be applied to any signal, because, as you already mentioned, it is unclear what it represents.


Need help deriving the energy of a signal I am solving old exam problems in preparation of my exam in signals, and I'm having trouble with a question.

For a continuous aperiodic signal with the spectrum $X(\omega)=\exp(-\omega^2)$
I am to derive the energy of the signal. I have made it this far, in agreement with the solution manual:

$W = \frac{1}{2\pi}\int_{-\infty}^{\infty}
\exp(-(\sqrt{2}\omega)^2)\,\mathrm{d}\omega$

The assignment then suggests substituting in the error function

$\mathrm{erf}(z) = \frac{1}{\sqrt{\pi}} \int_{-z}^z \exp(-y^2)\,\mathrm{d}y$

and exploiting that $\mathrm{erf}(z) \rightarrow 1$ towards infinity, which makes good sense, but the result in the solution manual before reducing the expression is

$W = \frac{1}{2\pi} \sqrt{\frac{\pi}{2}}
\mathrm{erf}(\sqrt{2}\omega)\Big\rvert_{\omega \rightarrow \infty}$

I don't understand where the $\frac{1}{2}$ in the square root comes from, so I would appreciate if anyone could help me understand this.

The factor just comes from the substitution $y=\sqrt{2}\omega$ in the integral defining the error function:

$$\begin{align}E_x=\frac{1}{2\pi}\int_{-\infty}^{\infty}|X(\omega)|^2d\omega&=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-2\omega^2}d\omega\\&=\frac{1}{2\pi}\frac{1}{\sqrt{2}}\int_{-\infty}^{\infty}e^{-(\sqrt{2}\omega)^2}d(\sqrt{2}\omega)\\&=\frac{1}{2\pi}\frac{1}{\sqrt{2}}\int_{-\infty}^{\infty}e^{-y^2}dy\\&=\frac{1}{2\pi}\sqrt{\frac{\pi}{2}}\lim_{x\to\infty}\textrm{erf}(x)\\&=\frac{1}{2\sqrt{2\pi}}\end{align}$$


How to detect whether the object has triangular shape? The task is to detect arrow from medical images. After segmenting, arrow-like objects can be connected with other connected components (noise). In this context, detecting a part of the object that characterizes an arrowhead would be a solution since template based shape descriptor does not perform well.
Please help!

enter image description here

I think you can find the orientation of the arrow using principal component analysis (or fitting an ellipse and taking major/minor axis directions).

To locate the arrow head, use distance transform. You'll have to preprocess the image, like flood-filling the region. Because the arrow head base is usually the widest part of the image, you should get a peak there.

Now you can use this point along with the two orthogonal vectors calculated earlier. In the image below I've translated those vectors to go through the distance-transform peak(I didn't do any calculations, the image is just for illustration).

The blue axis (second component in PCA or the minor axis of ellipse) now divides the image into two halves, and if you take projection of each half on to the blue axis, the half that contains the arrow head should give you a triangular profile (blue and orange profiles in the image below).

EDIT

Another approach would be to use the distance transform peak to isolate the neighborhood of the arrow head. If we threshold the distance transformed image by some fraction of the peak value, for the arrow head we should still get a triangular shape. Now we can prepare a mask centered on the peak to extract this region. Below is a simple Matlab code that does this for a specific image:

im = imread('arrow.jpg');
gray = rgb2gray(im);
bw = im2bw(gray, graythresh(gray));
di = bwdist(bw);
[mx, mxidx] = max(di(:));
[ypeak, xpeak] = ind2sub(size(di), mxidx);
bw2 = di > mx*.5;
se = strel('disk', double(ceil(mx*1.5)), 0);
mask = zeros(size(bw2));
mask(ypeak, xpeak) = 1;
mask = imdilate(mask, se);
figure,
subplot(1, 3, 1), imshow(bw), title('binary')
subplot(1, 3, 2), imshow(bw2.*di, []), title('thresholded dist')
subplot(1, 3, 3), imshow(mask.*di.*bw2, []), title('masked')

And the result is:


Confirmation on how to calculate phase differences at every frequency point between two wideband signals? I have a problem that I thought was going to be simple, but it has become surprisingly stubborn and I question my method...

I have used the method (described below), but I wanted to confirm that this was the right approach, and see if there was another method I could use.

What I basically have is two wideband signals, (on the order of 100 Mhz BW), in quadrature. So I have the I and Q of those two signals. What I am trying to find is the phase difference between those two signals as a function of frequency.

The method I am currently trying is the following:

I take my quadrature complex signals, and window each of them with a hamming window 3 times over.
I take the FFT of the windowed complex signals, at some arbitrary length, say 2048. Call the results $S_1(w)$ and $S_2(w)$.
I multiply $S_1(w)$ by the conjugate of $S_2(w)$.
I take the argument out of the complex exponential vector from the above step (inverse tangent of imag over real). This gives me the phase differences between the two signals as a function of frequency.

My questions on this method:

Is this the right approach for this simple type of problem?
When I tried this method for a simulation with 1, 2 or a couple tones, I can get good results. I find the index of where my tones exist, (look at the abs of the FFT), I take that index, and look at the $\Delta$phase value from step 4 above. This gives me the right phase.

However, when I look the $\Delta$phase values for frequency indicies where said those frequencies dont exist, I get wild answers. Why is this? I realize those frequencies dont exist but shouldn't the phases just be 0 - 0 = 0 in that case?

Based on this (with a couple tones), can I readily extend this to a wideband signal and also expect good results?
My last question regards the nature of the phase delay here - physically what is happening is that one of the signals (optical) is going through a device that delays different frequencies differently. Of course the delays are manifested as phase offsets between all frequencies... but at the same time there is a time delay associated with this as well. Even if we were to measure the $\Delta$phase between the two signals properly, wouldnt it still be impossible to measure the time delay between them at those frequencies? I can have a time delay of 10 wavelenghts, but a phase offset of 0.

Thoughts?

Thanks in advance!

EDIT FROM COMMENT FEEDBACKS:

I think I am clearer now as to what my question really is:  For a case where you have well separated tones and high SNR, my method can work.

However, what to do in the wideband case? In this sense, suppose $s_1[n]$ and $s_2[n]$ both contain two tones that are very close to each other in the frequency space. But your FFT can only give you so much frequency resolution. So instead of applying my method to two peaks, you only have one peak to go on. Now what?

What relation would the $\Delta$phase of that peak have to the actual $\Delta$phase of the  'true' peaks? Is it their average?
One method I am thinking of is correlating the two signals in the time domain, and working backwards to get the phase offsets per frequency taking advantage of the fourier translation property, although part of me thinks it might not work since there is no new information, just looking at the problem from a different domain.
Another thought I had for the wideband case would be to put the whole thing through many narrow band filters, and then apply my original method. Is there a name for something like this?

Thanks again

What you seem to want is the phase response of the device that your optical system is passing through. Based on that, your high-level approach makes intuitive sense overall; I'm not sure why you chose a Hamming window and applied it three times. In any event, though, if you have a linear discrete-time system with (complex-valued) frequency response $H(e^{j\omega})$, then:

$$
Y(e^{j\omega}) = H(e^{j\omega})X(e^{j\omega})
$$

where $X(e^{j\omega})$ and $Y(e^{j\omega})$ are the discrete-time Fourier transforms (DTFTs) of the input and output signals, respectively. If you express the above in polar format (magnitude and phase), then the phase response of the system is just:

$$
\angle H(e^{j\omega}) = \angle Y(e^{j\omega}) - \angle X(e^{j\omega})
$$

that is, the phase response of the system just expresses the amount of phase shift that the system applies to the input signal as a function of (normalized) frequency $\omega$. Your approach takes advantage of the fact that for two complex numbers $X(e^{j\omega})$ and $Y(e^{j\omega})$:

$$
\angle Y(e^{j\omega}) - \angle X(e^{j\omega}) = \angle \left(X^*(e^{j\omega})Y(e^{j\omega})\right)
$$

So, given the DTFTs of the system's input and output signals, you can straightforwardly calculate the phase response that you seek. The devil is in the details of how one would obtain the DTFT representations of these signals. The nice part about the DTFT is that it is continuous in frequency; the bad part is that the DTFT sum is infinite in length. So, if you really want to evaluate the phase response at arbitrarily many frequencies, you better be patient, because you'll need to collect an infinite amount of data first.

The discrete Fourier transform (DFT) is often used as an approximation to the DTFT, which is valid under certain conditions. If the signals $x[n]$ and $y[n]$ are finite in length (which for a practical scenario like yours, they are), then the DFT can be thought of as evaluation of the DTFT on the finitely-spaced frequency grid $\frac{2\pi k}{N}, k = 0, 1, \ldots , N-1$ (where $N$ is the signal length). So, the larger you make $N$, the smaller your sample spacing across the frequency axis.

Most physical devices will have a frequency response that is smooth to some degree (taken with the caveat that I'm not knowledgeable on optical device physics), so you should be able to find a DFT length that gives a decent approximation to the DTFTs for the input and output signals of your system. You can then apply your approach to estimate the phase response of the device.

You noted in your question that you see incorrect results for the estimated phase response for DFT bins in the output signal that didn't contain any content in the input signal; this is to be expected. Your assumption that you would just get zero phase for those bins is not likely to hold in a practical scenario; noise present in the sampled signals, spectral leakage, or other interference can cause you to observe some small amount of energy at frequencies that you wouldn't expect it. These extraneous sources can introduce distortion between the input and output that is not due to the device's response, and for complex numbers with small magnitude, small amounts of additive noise can drastically affect the number's resulting phase.

Therefore, you can only consider the phase measurement valid at frequencies where there is significant spectral content in the input. Stated differently, if you want to make a precise phase measurement, you want high signal-to-noise ratio in the complex value that you're using to estimate the phase.


Decimating Polyphase Filter in Simulink I have read many of the articles on the website regarding a polyphase decimating filter structure.  I'm sure we are all familiar with the typical image shown in Fred Harris Multirate book and many other papers written on this topic below.

Say we have some transfer function H(z).. well the polyphase implementation works by decimating the filter into M parallel filters where fs/M is the desired output rate.  The noble identities tell us we can pull the decimation through the filter stages... and hence the delays at the input are replaced by the commutator switch which effectively sends samples in at rate fs/M.  Let N be the total number of taps in H(z)... the convenience of this is that we only need to filter with N/M taps per cycle which greatly reduces resources... ok I've done my homework.

Now my question.  I have seen many implementations of this structure in code but I can't seem to find any implementations in simulink...specifically with what exactly is the commutator switch?

Here is an implementation that I have created in simulink... although it clearly doesn't work properly...Say I have some filter that was originally 120 taps and I now want an output rate that is 1/4 of the input rate... I decimate my filter into 4 sub filters and arrange into a parallel structure as shown below:

The switch on the input has a counter attached to it because I'm thinking that I can cycle between inputs at each clock edge to simulate the commutator...but then I need to demux the data into a branch that I want.. however this isn't working because the demux isn't controllable...

What am I doing wrong here?  What is the block that simulates the commutator like shown in many of the block diagrams in the multirate papers we all know?  I realize that it's a memoryless switch but which simulink block ?

Just to follow up on this, I believe I figured out one way to do what I was looking for shown below...if anyone has any other ideas feel free to share too


upconverting a kHz/MHz signal to GHz I want to know how in practice we can take a kHz/MHz signal to a GHz band?

Suppose we have an FSK modulated signal at around 300 kHz (or 20 MHz). I want to take it to 3.5 GHz, so about 3499700 kHz higher (or 3480 MHz higher).

Should I build an oscillator with a frequency of 3499700 kHz (or 3480 MHz) to use it as a mixer?
Is it possible in practice to have a precise 3499700 kHz oscillator?
Is there any IC useful or should it be by PLL?
How?
Is there any other and better solution?

I want to know how in practice we can take a kHz/MHz signal to a GHz band?

You guessed right: a mixer.

  Should I build an oscillator with a frequency of 3499700 kHz (or 3480 MHz) to use it as a mixer?

Assuming these frequencies are the difference between your source and your target frequencies, yes.

But you don't use an oscillator as a mixer, you use its output to feed it to a mixer.

  Is it possible in practice to have a precise 3499700 kHz oscillator?

Define "precise"!

Precise is a relative term. If you build a house, your standard of length needs to be accurate to say 0.2 % (so that if your house is 10m long, its length error will be no more than 2cm). If you build an aircraft engine, you might want to be better.

With oscillators it's usually way, way more precise than that, even on cheap standards.

Technologically, RF devices typically use synthesizers to generate local oscillator (LO) frequencies from a reference oscillator, which might be running, for example, at 10 MHz. A synthesizer typically contains all the parts necessary to derive that frequency – a control loop for a voltage-controlled oscillator, that oscillator, some sort of frequency multiplier, filtering to actually get the multiple you want, control and compensation circuitry…

The precision and stability (two different things!) of the resulting oscillator usually mainly depends on the precision and stability of the 10 MHz reference, and to a smaller extent (at least in the frequency regions you mention) on the noise and precision of the control loops.

  Is there any IC useful or should it be by PLL?

This question makes no sense. All the ICs that might make sense here contain at least one form or another of a PLL.

  How?

Good!

  Is there any other and better solution?

To what?

You seem to want a mixer + filters + a synthesizer. Mixers and filters can be bought as modules for moderate prices (e.g. from minicircuits), synthesizers often either take the shape of signal generators in laboratory settings, or as synthesizer ICs + support circuitry in actual devices.


Implementing Continuous Wavelet Transform I need to implement the discretized continuous wavelet transform from scratch. Could someone please point me to useful papers and references available online for this?

In 1D, some of the standard references are:

Continuous wavelet transform with arbitrary scales and $O({N})$ complexity, A. Muñoz and R. Ertl\'e and M. Unser, Signal Processing, 2002
A fast approximation to the continuous wavelet transform with applications, Berkner, K. and Wells, R. O., Jr., 1997, Proc. Asilomar
Fast Quasi-Continuous Wavelet Algorithms for Analysis and Synthesis of One-Dimensional Signals, Maes, S. H., 1997, SIAM J. Appl. Math.
Comparison of algorithms for the fast computation of the continuous wavelet transform, Vrhel, M. J. and Lee, C. and Unser, M. A., 1996, Proc. SPIE
Fast algorithms for discrete and continuous wavelet transforms, Rioul, O. and Duhamel, P., 1992, IEEE Trans. Inform. Theory

For a GPU implementation:

Effiziente Berechnung der FWT auf Grafikkarten, Richard Maltan, Bachelorarbeit, 2014


Coefficient meaning in DWT I understand that approximate and detail coefficient represent the different signal bands. But what do the values mean and how are they used?

DWTs turn a signal into a list of coefficients, whose combined sizes are equal to or greater than the signal length. Coefficients are grouped in packs called subbands. Each subband gathers coefficients resulting from some band-pass filtering, followed by a subsampling.
They are not exactly narrow-band signals, but projections onto embedded subspaces. The "narrowband" signal for each subband could be recovered by inverse transforming those coefficients only.


Discrete-time Fourier Transform of the unit step sequence $u[n]$ From text books we know that the DTFT of $u[n]$ is given by

$$U(\omega)=\pi\delta(\omega)+\frac{1}{1-e^{-j\omega}},\qquad -\pi\le\omega <\pi\tag{1}$$

However, I haven't seen a DSP textbook that at least pretends to give a more or less sound derivation of $(1)$.

Proakis [1] derives the right half of the right-hand side of $(1)$ by setting $z=e^{j\omega}$ in the $\mathcal{Z}$-transform of $u[n]$, and says that it is valid except for $\omega=2\pi k$ (which is of course correct). He then states that at the pole of the $\mathcal{Z}$-transform we have to add a delta impulse with an area of $\pi$, but that appears more like a recipe to me than anything else.

Oppenheim and Schafer [2] mention in this context

  Although it is not completely straightforward to show, this sequence can be represented by the following Fourier transform:

which is followed by a formula equivalent to $(1)$. Unfortunately, they didn't take the trouble to show us that "not completely straightforward" proof.

A book that I actually didn't know, but which I found when looking for a proof of $(1)$ is Introduction to Digital Signal Processing and Filter Design by B.A. Shenoi. On page 138 there's a "derivation" of $(1)$, but unfortunately it is wrong. I asked a "DSP-puzzle" question to have people show what is wrong with that proof.]

So my question is:

Can anybody provide a proof/derivation of $(1)$ that is sound or even rigorous while being accessible for mathematically inclined engineers? It doesn't matter if it's just copied from a book. I think it would be good to have it on this site anyway.

Note that even on math.SE almost nothing relevant is to be found: this question has no answers, and that one has two answers, one of which is wrong (identical to Shenoi's argument), and the other one uses the "accumulation property", which I would be happy with, but then one needs to proof that property, which puts you back to the start (because both proofs basically prove the same thing).

As a final note, I did come up with something like a proof (well, I'm an engineer), and I will also post it as an answer some days from now, but I would be happy to collect other published or unpublished proofs that are simple and elegant, and, most importantly, that are accessible for DSP engineers.

PS: I do not doubt the validity of $(1)$, I would just like to see one or several relatively straightforward proofs.

[1] Proakis, J.G. and D.G. Manolakis, Digital Signal Processing: Principles, Algorithms, and Applications, 3rd edition, Section 4.2.8

[2] Oppenheim, A.V. and R.W. Schafer, Discrete-Time Signal Processing, 2nd edition, p. 54.

Inspired by a comment by Marcus Müller, I'd like to show that $U(\omega)$ as given by Eq. $(1)$ satisfies the requirement

$$u[n]=u^2[n]\rightarrow U(\omega)=\frac{1}{2\pi}(U\star U)(\omega)$$

If $U(\omega)$ is the DTFT of $u[n]$, then

$$V(\omega)=\frac{1}{1-e^{-j\omega}}$$

must be the DTFT of

$$v[n]=\frac12\text{sign}[n]$$

(where we define $\text{sign}[0]=1$), because

$$V(\omega)=U(\omega)-\pi\delta(\omega)\Longleftrightarrow u[n]-\frac12=\frac12\text{sign}[n]$$

So we have

$$\frac{1}{2\pi}(V\star V)(\omega)\Longleftrightarrow \left(\frac12\text{sign}[n]\right)^2=\frac14$$

from which it follows that

$$\frac{1}{2\pi}(V\star V)(\omega)=\text{DTFT}\left\{\frac14\right\}=\frac{\pi}{2}\delta(\omega)$$

With this we get

$$\begin{align}\frac{1}{2\pi}(U\star U)(\omega)&=\frac{1}{2\pi}\left[(\pi\delta(\omega)+V(\omega))\star (\pi\delta(\omega)+V(\omega))\right]\\&=\frac{1}{2\pi}\left[\pi^2\delta(\omega)+2\pi V(\omega)+(V\star V)(\omega)\right]\\&=\frac{\pi}{2}\delta(\omega)+V(\omega)+\frac{\pi}{2}\delta(\omega)\\&=U(\omega)\qquad\text{q.e.d.}\end{align}$$

Cedron Dawg posted an interesting initial point in this answer. It begins with these steps:

$$ \begin{align}
U(\omega) &= \sum\limits_{n=0}^{+\infty}  e^{-j \omega n} \\
&= \lim_{ N \to \infty } \sum\limits_{n=0}^{N-1}  e^{-j \omega n}\\
&= \lim_{ N \to \infty } \left[  \frac{  1 -  e^{-j \omega N}  }{ 1 -  e^{-j \omega }   } \right] \\
&= \frac{ 1 }{ 1 -  e^{-j \omega } } - \lim_{ N \to \infty } \left[  \frac{  e^{-j \omega N}  }{ 1 -  e^{-j \omega } } \right]
\end{align} $$

It turns out the term inside the limit can be expanded as follows:

$$\begin{aligned}
\frac{  e^{-j \omega N}  }{ 1 -  e^{-j \omega } }
={}
&\frac{1}{\sin^2(\omega)+(1-\cos(\omega))^2} \cdot \\
&\left[-\cos(\omega)\cos(N\omega)+\cos(N\omega)-\sin(\omega)\sin(N\omega)+ \\
j(-\sin(\omega)\cos(N\omega)+\cos(\omega)\sin(\omega)-\sin(N\omega))\right]
\end{aligned}
$$

The common factor outside the brackets can be expressed as:

$$\frac{1}{\sin^2(\omega)+(1-\cos(\omega))^2} = \frac{1}{4\sin^2(\omega/2)}$$

The real part inside the brackets also equals:

$$-\cos(\omega)\cos(N\omega)+\cos(N\omega)-\sin(\omega)\sin(N\omega)= 2\sin(\omega/2)\sin[\omega(-N+1/2)]$$

On the other hand, the imaginary part can be rewritten as:

$$-\sin(\omega)\cos(N\omega)+\cos(\omega)\sin(\omega)-\sin(N\omega)=-2\sin(\omega/2)\cos[\omega(-N+1/2)]$$

Rewritting the original term we get that:

$$\begin{align}
\frac{  e^{-j \omega N}  }{ 1 -  e^{-j \omega } }
&=\frac{2\sin\left(\frac \omega 2\right)}{4\sin^2\left(\frac \omega 2\right) }
\left( \sin[\omega(-N+1/2)] - j\cos[\omega(-N+1/2)]\right) \\
&=-\frac{\sin[\omega(M+1/2)]}{2\sin\left(\frac \omega 2\right)} - j\frac{\cos[\omega(M+1/2)]}{2\sin\left(\frac \omega 2\right)}
\end{align}
$$

where I used $M=N-1$ and the limit stays unaffected as $M\to \infty$ as well.

According to the 7th definition in this site:

$$\lim_{M\to \infty} -\frac{1}{2\sin(\omega/2)}\sin[\omega(M+1/2)] = -\pi \delta(\omega)$$

So far we have that:

$$\lim_{ M \to \infty }   \frac{  e^{-j \omega (M+1)}  }{ 1 -  e^{-j \omega } } =-\pi\delta(\omega)
-j\lim_{M\to \infty} \frac{\cos[\omega(M+1/2)]}{2\sin(\omega/2)}
$$

If we could prove that the second term on the right of the equality is $0$ in some sense, then we are done. I asked it at math.SE and, indeed, that sequence of functions tends to the zero distribution. So, we have that:

$$ \begin{align}
U(\omega) &= \frac{ 1 }{ 1 -  e^{-j \omega } } - \lim_{ N \to \infty } \left[  \frac{  e^{-j \omega N}  }{ 1 -  e^{-j \omega } } \right] \\
&=\frac{ 1 }{ 1 -  e^{-j \omega } }+\pi\delta(\omega)
+j\lim_{M\to \infty} \frac{\cos[\omega(M+1/2)]}{2\sin(\omega/2)} \\
& =\frac{ 1 }{ 1 -  e^{-j \omega } }+\pi\delta(\omega)
\end{align} $$


Deconvolve a FIR-filter using limits? Or: Are limits distributive with the inverse fourier transform? Given $y[n] = f[n] * h[n]$ where $y[n]$ and $f[n]$ are two one-dimensional discrete  signals that are given, find out:

if $h[n]$ is a FIR-filter
if it is a FIR-filter, what its kernel is.

Obviously, this can be done by doing deconvolution: $h[n] = y[n] *^{-1} f[n]$.
Performing deconvolution (written here as $*^{-1}$) in the time domain is of course computationally expensive. However, it is possible to go to the frequency-domain instead:

$$h[n] = y[n] / f[n] = \mathcal{F^{-1}}(Y ./ F) $$
(Where $ ./$ is pairwise division and $Y$ and $F$ are the (discrete) fourier transformations of $y[n]$ and $f[n]$ respectively, that is $Y = \mathcal{F}(y[n]$), $F = \mathcal{F}(f[n])$).

Many existing algorithms use this property, since pairwise division has a lower time-complexity) than direct deconvolution.

However, these algorithms claim to struggle with zero-valued frequencies: If $F$ contains a zero, then we will divide by zero during the pairwise division.

My question is: Is it possible to calculate the outcomes using limits here? After all:

\begin{align}
\lim_{k\to 0^+} \frac{x}{k} &= \lim_{l\to \infty} l \tag{$0^+$}\\
\lim_{k\to 0^-} \frac{x}{k} &= \lim_{l\to -\infty} l  \tag{$0^-$}\\
\mathcal{F}(x) &= \int_{-\infty}^{\infty} f(x)\ e^{-2\pi i x}\,dx\tag{FT}\\
\mathcal{F^{-1}}(x) &= \int_{-\infty}^{\infty} f(x)\ e^{2\pi i x}\,dx\tag{IFT}\\
\int_{-\infty}^{\infty} f(x)\,dx &= \lim_{h \to 0} \sum{\frac{f(x)+f(x-h)}{2}} \tag{sum limit integral}\\
\end{align}

So the fourier transform (and the inverse fourier transform) are defined in terms of integrals, which are themselves defined in terms of limits.

Limits distribute over summation, so it seems to me that it would be possible to simply use IEEE 754 floating-point numbers (which follow the limit-rules for division by zero) to find a correct answer, which would allow for a very fast computer-implementation of a deconvolution algorithm.

As an aside, I wonder if the fact that $h[n]$ is a FIR-filter (and linear, time-invariant and causal) is useful in some way.

Is my reasoning correct, or is it flawed somewhere?

Integral definition

So the fourier transform (and the inverse fourier transform) are defined in terms of integrals, which are themselves defined in terms of limits.

Nope, the integrals are not just the limits of sums here. That would be true for Riemann integrals over smooth functions (try doing that difference quotient at a place where $f$ isn't continuous). It's not how the Lebesgue integral that we need to use when defining the Fourier transformation is defined.
And, especially, if you understand the DFT as a "special use case" of the continuous FT, then you have to assume periodic repetition of the signal vector you're observing.
Periodicity leads to line spectra. These are, by definition, not continuous.
So, no, your last equation isn't right.
Fourier Transform definitions
Your equations $(\text{FT})$ and $(\text{IFT})$ aren't how we write "The Fourier transformation of $f(x)$ is ...": Your left side of the equation doesn't say you're transforming $f(x)$; it doesn't mention $f$ at all!
One common way I see the "Fourier transform of $f$" (which is a function of the frequency $\omega$) being written is
\begin{align}
\mathcal{F}\{f\}(\omega)
\end{align}
and while we're at it, the right side of these equations is wrong, too: you're missing the $\omega$ (or whatever target-domain free variable you want to use) in the exponent.
Limit $l$
Now, the first two equations:
\begin{align}
\lim_{k\to 0^+} \frac{x}{k} &= \lim_{l\to \infty} l \tag{$0^+$}\\
\lim_{k\to 0^-} \frac{x}{k} &= \lim_{l\to -\infty} l  \tag{$0^-$}
\end{align}
stand a bit on their own here, because you don't explain what $l$ is meant to be, but:
The right hand side of $(\text{$0^+$})$ and $(\text{$0^-$})$ don't exist. I mean, you literally write "when we not fix $l$ to any actual value, but let it simple pine for infinity, what's the actual value?". So, that actual value then doesn't exist.
Finding a substitute value of a function in a singularity

Is my reasoning correct, or is it flawed somewhere?

Yeah, at least the way it's written down, it's mathematically flawed.
Of course, you're right from a functional analysis point of view, your quotient function simply has a couple of singularities. That's not a problem per se!
Some of these singularities can actually be solved by finding a series on the environment of that point that is equal to the quotient in all but that point, and converges, and then "fix" the value at that point. However, that puts a huge restriction on your quotient function: It needs to be holomorphic in the neighborhood of that singularity for that series to exist. Sadly, quotients of holomorphic functions are only holomorphic where the denominator isn't zero. So, bad luck going the limit route!
What you're looking for is the residue theorem! I'll really quickly outline the idea:

your quotient function is generally holomorphic, aside from a discrete set of points (where your denominator is 0 or your numerator goes towards infinity).
If a function is holomorphic, than it has the (funky) property that if you have a "disk" (a non-intertwining surface with a closed boundary)  on which it is defined, then you can fully infer the function from the values on the boundary of that disk. You might have met that, if you're an EE, used in Stoke's theorem (in electrodynamics)
So, you just define a small disk around your singularities with radius $\gamma$, and let that radius get smaller. That should define how your quotient should behave in the singularity.

Problems
The above thing is a very analytical way of describing a system. Yes, if you have an analytic description of your quotient (hint: you do, see what sampling does and how we mathematically define the DFT as finite sum over holomorphic functions), you can apply that knowledge.
The practical problem arises from the fact that your quotient is based on an observation. Observations are noisy. Where an observed quality is close to 0, the noise component becomes dominant. If you divide by that, you get a huge amplification of noise!
So, whilst nice analytically, that's not the way you want to go in the real world. That's why all the deconvolution algorithms (there's actually quite a few) apply assumptions on the system being convolved with, and minimize some error term or maximize some likelihood term.


Summation of two IIR filters Actually my question is how to get the weighted average of two or three IIR filters given their transfer functions, which is nearly equal to this question. How to calculate the total transfer function of the summation of two IIR filters? For example, I have two filters as following
$$H_1(z) = \frac{b_0+b_1z^{-1}+b_2z^{-2}}{1+a_1z^{-1}+a_2z^{-2}}$$
$$H_2(z) = \frac{b_0'+b_1'z^{-1}+b_2'z^{-2}}{1+a_1'z^{-1}+a_2'z^{-2}}$$
I can directly solve the sum of these two fractions as $H(z) = H_1(z)+H_2(z)$, but the problem is that my filter orders are generally around dozens, I have no idea how to do it.
I also know that I can put a signal through this two filters and sum the output together, but it requires double computational complexity.
Is there any other way to solve this problem? Thanks!

In theory you just do it the same way. The transfer function of a weighted sum of filters with weights $w_i$ is just
$$H_{total}(z)=\sum_iw_iH_i(z)\tag{1}$$
where $H_i(z)$ are the individual transfer functions.
However, in general the order of the total transfer function is the sum of the orders of the individual filters, unless some filters share the same poles. Consequently, the resulting filter order can become too large for a practical implementation. If the orders of the individual filters "are generally around dozens" then you have a problem, even using floating point arithmetic.


Understanding the phase in DFT Although using regularly the FFT algorithm to compute DFTs, I don't really understand how the phase part works.

e.g. in Matlab:

Fs=1e3;
t=linspace(0,1,Fs);
f=1;
x=sin(2*pi*f*t);
X=fft(x);
Xphase=angle(X);
Xmag=abs(X);

I know that "phase" can mean a shift (an angular) shift of a wave with respect to another. But I can't really say what the values in each point of Xphase mean. Can someone explain ? e.g. I was trying with this simple example of a pure sine wave (which shouldn't be dephased I guess?)

Xphase contains values between -3.1385 to 3.1416 so roughly from -180 degrees to 180 degrees... why?

Suppose you have a sinusoidal that has a whole number of cycles ($k$) in your DFT frame containing $N$ sample points.  It can be parameterized like this:

$$ x[n] = A \cos \left( \left( k\frac{2\pi}{N}\right)n + \phi \right) $$

If you take the $1/N$ normalized DFT of this (FFT is a DFT that is computed efficiently), all the bins will be zero except for bins $k$, and $(N-k)$.  With MATLAB, bin $k$ occurs at index $k+1$.

$$ X[k] = \frac{A}{2} e^{i\phi} $$

and

$$ X[N-k] = \frac{A}{2} e^{-i\phi} $$

So, you can see, in the ideal case of a pure tone with a whole number of cycles in the frame, the phase angle of the DFT bin corresponds directly to the phase argument in the signal.

The values from $-\pi$ to $\pi$ are by convention and are measured in radians.  This range covers every possible angle.

If you don't have a whole number of cycles, you can find my simplified bin value formulas here:  https://www.dsprelated.com/showarticle/771.php

https://gizmodo.com/pentagon-ordered-to-tell-congress-if-it-weaponized-tick-1836391549

In the time domain, a $2\pi$ change in the phase represents the shift of an entire cycle, which can also be considered one wavelength.  Shift it by $\pi$ (popularly known as a 180 degree shift), and you effectively invert it.  Shift it back and forth by $\pi/2$ and sine becomes cosine and vice versa.

Maybe this will spark an understanding:

$$ cos( ft + \phi ) = \cos( ft ) \cos( \phi ) - \sin( ft ) \sin( \phi ) = a \cos(ft) + b \sin(ft) $$

So, fiddling with the phase adjusts how much cosine vs the sine is in the tone within that reference frame.

What I described above is the bedrock connection between the phase in the time domain and the phase in a DFT bin for real valued signals.

The relationship between this shifting in the time domain and the angle in the corresponding bin is one-to-one for sinusoidals with a whole number of cycles in a DFT frame.

$$ a = \cos(\phi) $$

$$ b = -\sin(\phi) $$

$$ \frac{b}{a} = -\frac{\sin(\phi)}{\cos(\phi)} = -\tan(\phi) $$

Mix in an $i^2$ and you got bin interpretation.


Why histogram equalization fail in enhance images for bimodal histogram? I'm on histogram equalization process. I have quite well understood the math but i don't understand why an histogram on bimodal image fails in image enhancement and what techniques can be used in this case (some kind of local/adaptive histogram?).

Yes it is true. A simple histogram equalization may fail for certain image intensity distributions. In your case there are some bins with excessive number of pixels. For other methods, try for example gamma and contrast corrections.


Super resolution methods I am working on the project where my goal is to create SR image from 25 LR images that are precisely shifted by 0.2px (both in horizontal and vertical direction ). I also know that the captured object is a back and white grid ( something like chessboard)
I am looking for some precise SR method that will take advantage of the precise shift of the images. The emphasis is on precision not on the speed of the method. I have read a lot of the papers on that topic and found out that there is a lot of methods and it is difficult for me to decide which method suits my task.

If you have experience with SR please suggest some methods that find good for this task.

You'll need to describe the density ("kernel") with which light from different angles contributes to each (physical) pixel's amplitude. Expect something rather gaussian.

Now, shift that kernel by 0.2 pixel width, repeat. You'll find a Max-Likelihod estimator for the high-density picture simply by finding a system of equations that describe how much the values of the 5 pixel positions whose densities overlap that position contribute.

If you want a very intuitive first approach, look for the kdeplot tool of the python seaborn library. The website has a nice gallery.


Inverse Fourier Transform of piecewise function Objective : Compute $y(t)$ from $Y(j\omega)=H(j\omega)X(j\omega)$ where :
$$
x(t)=\left(\frac{\sin(2t)}{\pi t}\right)^{2}
$$
and
$$
H(j\omega)=\begin{cases}e^{-j\omega}&\text{if $|\omega|<4$}\\0&\text{if otherwise}\end{cases}
$$
what I did is find FT of $\left(\frac{\sin(2t)}{\pi t}\right)\left(\frac{\sin(2t)}{\pi t}\right)$ as follow :
$$
\mathcal{F}\{s\cdot s\}(t)=\frac{1}{2\pi}[S*S](j\omega)
$$
We have : $$S(j\omega)=\mathcal{F}\left\{\frac{\sin(2t)}{\pi t}\right\}=\begin{cases}1&\text{if $|\omega|<2$}\\0&\text{if otherwise}\end{cases}$$
I have computed the modulation :
$$
X(j\omega)=\frac{1}{2\pi}[S(j\omega)*S(j\omega)]=\frac{1}{2\pi}\int_{-\infty}^{\infty}S(j\Omega)S(j(\omega-\Omega))\;\text{d}\Omega=\begin{cases}\frac{\omega}{2\pi}&\text{if $0<\omega<2$}\\\frac{2-\omega}{2\pi}&\text{if $2<\omega<4$}\\0&\text{if otherwise}\end{cases}
$$

The first issue I have is multiplying $H(j\omega)$ with $X(j\omega)$ since I have two piecewise functions. The second issue is how to thus compute its inverse Fourier transform.

I would hope someone can please help me. Thank you

So you know that in the frequency domain the function
$$s(t)=\frac{\sin(2t)}{\pi t}\tag{1}$$
is a rectangle with height $1$ for $|\omega|<2$. Now think graphically about convolving $S(j\omega)$ with itself. Can your result be right? You found that $X(j\omega)$ is only non-zero for positive frequencies. That would correspond to a complex-valued time domain function. So check that result again. The result of the convolution of two rectangles is basic and you should be able to see the result graphically.
After having figured out the convolution of $S(j\omega)$ with itself, it should be easy to obtain the bandwidth, i.e. the frequency $\omega_c$ above which the spectrum of the result is zero. With that knowledge it's straightforward to apply $H(j\omega)$ if you know the time-shifting property of the Fourier transform.


Regression vector size for prediction, reconstruction and filtration with adaptive filters I am working with adaptive filters and similar adaptive models (mainly with gradient adaptation) for a few years. I and my colleagues always struggle to find out the correct size of regression vector.

So far, as I find out by the hard way (experience):

for filtration mostly works well (low MSE) with something about 10-50 last samples, this number changes according to data and used step-size / algorithm.
for prediction the optimal size also depends on data and learning rate. Longer windows (regression vectors) commonly works well and it can be beneficial to skip samples, example:

Example ad 2.

For data where is one really dominant frequency (but changing over time) works for me the vector

$[u(k), u(k-15), u(k-30), u(k-45), u(k-60)]$ (size of 4 samples)

with the same MSE as

$[u(k), u(k-1), ..., u(k-60)]$ (size of 60 samples)

it seems that only what matters is to cover the size of data about 2 last waves of with the most dominant frequency, no matter how much samples was used.

My question

What is applicable theory for this problem (at least possibly)? I am pretty sure there must be some data attributes what can tell how big regression vector is important to get the most from the data. But I do not really now what should I study, or how to search for it.

Any hint or research direction recomendation is welcome. Thanks in advance.

Suppose that you are adapting $w$ to minimize $\text{E}(y[n]-w[n]*u[n])^2$ where
$$y[n]=h[n]*u[n]+\nu[n]$$
$y[n]$ and $u[n]$ are known and $\nu[n]$ is an additive noise component.
With a long enough FIR filter you can model any linear system $h$ with any accuracy. You can opt for a shorter filter length in expense of limited modelling accuracy. A key parameter is the sampling rate of the signals. If the signals are over sampled unnecessarily then the required length for $w$ increases.

Suppose that we settle for a certain modelling accuracy and we want to find $w$. If all the samples of $u[n]$ are available and $h$ remains constant you can find an estimate of $h$ namely $w$ using a least square estimation. This estimation will have a limited accuracy. But it is the best you can do. This error is referred to as the Least square error which here I (a bit wrongly) consider the same as Mean Square Error (MSE error).

Since $h$ may change, we need to continuously estimate $h$, and that is why we are using the adaptive filters. Let's say that we use the LMS adaptation to find $w$.
Now on top of the MSE error we have extra error owing to the way that the LMS adaptation is performed. This error is called  misadjustmentis and it is caused by the limited window where $u[n]$ is used in the adaptation. If the LMS step size $\mu$ approaches zero we are including all sample of $u[n]$ in estimation, but then there is no adaptation. It can be proved (Take a look at any adaptive filters book) that misadjustment is proportional to

The filter length
The step size

The other issue is the tracking behavior of the filter. Since $h$ is changing we like to adapt the filter as fast as possible. To increase the adaption speed we can increase $\mu$. To keep the filter stable we should have
$\mu<\frac{1}{\lambda_{\text{max}}}$
where $\lambda_i$ are the eigenvalues of the covariance matrix of $u[n]$.
The situation is worse if the sample of $u[n]$ are correlated (apparantly your scenario). In this case the covariance matrix has different eigen values and the adaptive filter converges with different modes where the time constant $\tau_i$ for each mode will be
$$\tau_i=\frac{1}{2 \mu \lambda_i}>\frac{\lambda_{\text{max}}}{2  \lambda_i}$$
The ratio of largest to smallest eigenvalue $\frac{\lambda_{\text{max}}}{\lambda_{\text{min}}}$ increases as the samples of $u[n]$ become more correlated and the LMS filter becomes slower. Unfortunately the speed cannot be increased idefinitely by increasing $\mu$. Since $\mu$ is limited from above by $\frac{1}{\lambda_{\text{max}}}$.

Back to your case:
It seems that your data is highly correlated that is why removing 14/15 of samples does not change the performance that much.
Throwing away samples of $u[n]$ can help in your case for the following reasons:

You are decreasing the length of the part that you are adapting. This can decrease the misadjustment.
You are whitening your data by down sampling it. This decreases the eigen ratio and improves the tracking behavior of the filter.
Another issue with highly correlated $u[n]$ is the potential divergence of filter at frequencies where $u[n]$ has almost no component. If $u(e^{j\omega_i})=0$ then $h(e^{j\omega_i})$ can go to infinity if the filter length is infinite. The shorter the filter the least likely that this occurs. This problem however can be solved by introducing leakage.

The negative impact of throwing way samples of $u[n]$ is that your are imposing  a certain assumption on $w$. By using only $[u(k), u(k-15), u(k-30), u(k-45), u(k-60)]$  you are assuming that $w$ takes a form like:
$$w=[w_0~...~0~w_{15}~0~...0~w_{30}~...~w_{60}]$$
This limits the ability of $w$ to model an arbitrary $h$.


How to construct an unknown signal from sparse periodic samples? I have an unknown black-box that is generating a deterministic continuous periodic signal and the period is variable between 1 to 30 seconds.

When the black box is turned on, the period is fixed, but unknown.
Sampling rate and signal period are not synched, and it is highly unlikely to be integer coefficient of one another.

All that is available to me is sparse (contains missed) samples with period of 100 ms to 10 ms.  (The parameter of source is fixed every time, for example I should construct a 3s signal from periodic 10ms samples). There is no limit on the duration of sample gathering but samples may be missed. Signal is smooth and does not change fast and slowest sampling rate satisfies nyquist rate.

Is there a method to construct a single period of signal and enhance its resolution as time goes by?

Nonuniform Sampling

Assuming you know what the actual sampling rate is, and exactly which samples are missed, it is possible to reconstruct the signal over a time interval T as long as the bandwidth W of the signal still satisfies the Nyquist criterion: W < N/T/2, where N is the number of samples you do have.

Visualize it this way:
A signal of bandwidth W over a time interval T can be decomposed into a
linear combination of N sinc functions (sinc is shorthand for sin x/x),
where N = 2×W×T.
The sinc functions are equally spaced over the interval in question and
their widths are all the same.
At the instant where each sinc pulse reaches its peak, all of the other
sinc functions have a value of zero.

If you sample the signal at the instants corresponding to the peaks of
the sinc functions (tA0, tA1 and tA2),
each sample directly yields the amplitude of the corresponding sinc pulse
(A0, A1 and A2).
The original signal can be reconstructed by creating new sinc pulses of the
correct amplitudes and overlaying them, usually by using a digital or
analog filter with the correct brick-wall frequency response.
The reconstructed signal will be described by the equation:

Signal(t) = A0 sinc (t – tA0)
     + A1 sinc (t – tA1)
     + A2 sinc (t – tA2)

If you take N samples of the signal at other instants, each sample will be
a linear combination of all of the sinc functions.
Each sinc pulse will contribute an amount to the sample based on the
time interval between the sample and the center of the sinc pulse.
It is still possible to calculate the amplitudes of the sinc pulses, but
now it is necessary to solve a solve a system of linear equations in order
to do so.

If you have samples B1, B2 and B3 taken at times tB1,
tB2 and tB3 as shown above, you can write the
following equations:

B1 = A0 sinc (tB1 – tA0)
     + A1 sinc (tB1 – tA1)
     + A2 sinc (tB1 – tA2)

B2 = A0 sinc (tB2 – tA0)
     + A1 sinc (tB2 – tA1)
     + A2 sinc (tB2 – tA2)

B3 = A0 sinc (tB3 – tA0)
     + A1 sinc (tB3 – tA1)
     + A2 sinc (tB3 – tA2)

Given that you know B1, B2 and B3, and all of the time values, it's a
straightforward matter to solve this system for A0, A1 and A2.
Once you have these coefficients, you can plug them into the reconstruction
equation given previously to reconstruct the original signal.

Note that in the above explanation, I have implicitly assumed that the
uniform sampling interval (tA1 – tA0) equals
π, in order to simplify the diagrams and
equations.
The same argument holds for other sampling intervals (and bandwidths) as
well.

Averaging Multiple Periods

Using multiple periods of the signal to improve the reconstruction requires that you can precisely correlate those periods in the time tomain. Once you have done the reconstruction described above, it should be possible to precisely locate features such as peaks and zero-crossings, even if these features fall between samples. Use this information to time-align multiple periods so that they can be averaged together to reduce errors such as noise, quantization error or timing jitter.


Methods to evaluate the quality of an IIR Filter fixed-point implementation I'm implementing a scattered look-ahead IIR filter in Matlab/VHDL, basically I convert an order-2 IIR filter into an order-4 IIR in order to take advantage of the pipelining of an FPGA. (reference http://people.ece.umn.edu/users/parhi/SLIDES/chap10.pdf)

Anyway, I'm trying to evaluate how many bits I need in every stage and I was wondering if there is a metric to evaluate the "quality" of the quantization.

For example an error of 0.1% in the transient of a step response doesn't matter much but a a neverending limit cycle of 0.1% after the transient would be a problem.

I was thinking of using a metric like m = error^2 * n

where error is the difference between the reference double implementation and my fixed-point implementation and "n" is the sample number so an error at sample 100 would weigh more than error at sample 1.

I would need to balance quantization errors with FPGA resources and FPGA speed of course.

Are there other metrics? Am I on the right track?

Any type of test result will be extremely dependent on the implementation, the filter coefficients and the actual input signal. Typically you would start with a detailed mathematical analysis of your filter implementation and calculate the transfer function from the input to each of the state variables. Based on either sine-wave scaling or "absolute sum of impulse response scaling" you can determine the max value of your states and determine the scaling that prevents clipping. The number of bits is then determined by your required SNR and how the noise propagates from the state variables to the output.

Next your check your rounding behavior: Safest is "round towards zero". That leaves a little of of SNR on the table, but it prevents any type of limit cycles.

Finally you test it. You either construct the worst case signals and or create a set of test signals that covers anything that can happen in your application including some of the corner cases. Than, indeed, a simple "RMS error" metric is a good choice.


How do real-time convolution plugins process audio so quickly Okay, so in Logic Pro I can load up a Space Designer plugin (convolution reverb) with an impulse that's 9.1 seconds long, turn my mic on, and get real-time convolution reverb as the mic records incoming audio.

My question is: How could this be possible? My I/O buffer size is 32 samples. Is my computer capable of performing the convolution of a 401,310 sample impulse with a 32 sample block of audio in under 0.0007 seconds? If so, what kind of black magic do I need to get my own convolution reverb to do this?

I should mention that I am trying to develop real time convolution in C with portaudio. It works great when it's not real-time, but it's way too slow in real-time (the convolution I mentioned above takes 0.8 seconds using frequency-domain convolution which is supposed to be faster for impulses of that size, no?

Anyway any insight would be helpful, as I would really like to understand how this kind of real-time processing is possible.

Real-time low-latency partitioned convolution reverb with a long impulse response works by dividing the impulse response into unequally sized partitions. The shortest partitions (blocks) are at the beginning of the impulse response, and the partition length grows towards the end of the impulse response:

Each partition length can be processed separately, giving for each a latency that is for fast Fourier transform (FFT) -based convolution perhaps a couple of times the partition length. Convolution outputs for all partition lengths are summed forming the final composite output. When partitioning the impulse response, the partition of a particular length must not appear before the time that equals the latency of convolution of that partition length. A time-domain finite impulse response (FIR) filter can be used for the beginning of the impulse response because that may be more efficient than doing the low-length convolution by FFT.

What is the optimal partitioning for a given impulse response length and the optimal split point between time-domain FIR and FFT convolution should be tested on the target platform. A challenging aspect of non-uniform partitioned convolution is scheduling of the work.

Non-uniform partitioned convolution was popularized by William G. Gardner, Efficient Convolution without Input-Output Delay, J. Audio Eng. Soc., Vol. 43, No. 3, 1995 March, available online at http://www.cs.ust.hk/mjg_lib/bibs/DPSu/DPSu.Files/Ga95.PDF.


What is meant by 20 dB signal-to-noise ratio? I'm reading a paper in which a discrete signal

$$x(n) = s(n) + w(n)$$

is considered. $s(n)$ is a known, deterministic series, and $w(n)$ is white noise with zero mean. The authors write that

  signals were generated with an SNR of 20 dB

What does this mean?

What does signal energy mean? There seems to be several ways of defining this, but no attempt to do so is done in this paper.
What does 20 dB SNR mean? I suppose that either $20=10\log{E_s/E_w}$, or $20=20\log{E_s/E_w}$, but am not familiar enough with these concepts to be sure.

A discrete signal is often interpreted as an analog voltage signal. The signal x(n) should be considered as a sampled version of x(t) with a unknown sampling frequency.

Q1:

So when we're talking about the energy of the signal, it is presumed that the voltage signal is obtained in a circuit with an 1Ω resistor in series, so that v(t)^2/R = P.

Q2:

Signal-to-noise ratio is used in many areas and is defined as the logarithmic power ratio of the signal and noise, i.e. SNR = 10*log(Ps/Pw). The other SNR = 20*log(Ps/Pw) is not true, because the Ps and Pw are power values. The latter is used when the signals are represented as "voltage" signals, i.e. SNR = 20*log(s(n)/w(n)) = 10*log(s(n)^2/w(n)^2).


Common Filter Types for Audio Applications I'm beginning to write a basic "virtual analog" synthesizer and want to implement a lowpass with variable cutoff frequency and resonance. Does anyone have insight into filter types commonly used to mimic the sound of traditional analog gear?

More generally, what are common filters types found in audio applications such as synthesizers, equalizers, etc. I took a basic DSP class many years ago but never gained an appreciation for design choices to consider for audio. I would appreciate any references along these lines.

a lowpass with variable cutoff frequency and resonance.

RBJ has a cookbook on how to do this:

http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

Here the cutoff frequency is called w0, and the resonance is called Q:

LPF:        H(s) = 1 / (s^2 + s/Q + 1)

            b0 =  (1 - cos(w0))/2
            b1 =   1 - cos(w0)
            b2 =  (1 - cos(w0))/2
            a0 =   1 + alpha
            a1 =  -2*cos(w0)
            a2 =   1 - alpha

If you want to vary these parameters continuously while passing signal, you are not supposed to use biquads (because updating the coefficients causes glitches?) and you should use state-variable filters instead:

http://www.earlevel.com/main/2003/03/02/the-digital-state-variable-filter/

This is also a 2nd-order lowpass with variable Q.

  Does anyone have insight into filter types commonly used to mimic the sound of traditional analog gear?

Generally any filter transformed to a digital IIR using the bilinear transform is an "analog-like" filter.  Bessel is an exception because the important property of a Bessel filter is its group delay, which is not preserved by the bilinear transform.  Bilinear transform attempts to preserve the frequency response primarily.


BIBO stability of $y(t)=\int_{-\infty}^{t}{x(\tau)d\tau}$ How can I prove that the LTI system with (output $y(t)$, input $x(t)$)

$$y(t)=\int_{-\infty}^{t}{x(\tau)d\tau}$$

is BIBO (bounded-input/bounded output) stable?

Since this is most likely homework, here is a hint.

Write the integral you have displayed in the form $\int_{-\infty}^\infty x(\tau)h(t-\tau) d\tau$ where you get to choose what the function $h(\cdot)$ is to make it all work out. Then, $h(t)$ is the impulse response of the LTI system, Do you know the criterion for BIBO stability of an LTI system in terms of $h(t)$? Does your $h(t)$ satisfy the criterion?

Edit in response to OP's comment and definition of BIBO stability:

Your understanding of the definition of BIBO stability is faulty.

  A system with input $x(t)$ and output $y(t)$ is said to be bounded-input bounded-output (BIBO) stable if it has the property that whenever $x(t)$ is bounded (that is, there exists a (finite) real number $M$ such that $|x(t)| \leq M$ for all real numbers $t$), then the output is also
  bounded (that is, there exists a (finite) real number $N$ (possibly dependent on $M$) such that $|y(t)| \leq N$ for all real numbers $t$).

Note that while $N$ might depend on $M$, it does not depend on $t$, the same $N$ must be an upper bound on $|y(t)|$ for all $t$.

In your case, with $y(t) = \int_{-\infty}^t x(\tau) d\tau$, what happens if $x(t) = u(t)$ for all real numbers $t$ (step input)? This is a bounded input. Does there exist a finite number $N$ such that $|y(t)| \leq N$ for all real numbers $t$?


How to add a particular noise on a signal? I have a simple ADC model (an input signal == sampled signal ). The adc should be modelled as with following values:

1.5 V (peak to peak).. meaning --> +1.5 V --> 0 dBFS ... how to do this ?
The noise on the signal should be -120 dBFS/Hz in bandwidth of 2 MHz.. .how to add this exact noise ?
The NFFT= 2^11 and fs= 200 MHz.
The input will be 10 MHz single tone, and then the FFT will be perfomed.

fs = 200e6;

% Adding Noise
n_en=1; % Noise On/Off (1/0)
BW= 2e6;
noise= -120; %dBFS/Hz
SNR=abs(noise + 10*log10(BW));

N= 2^11;
Df_fft = fs/N;
t = (0:N-1)/fs; t = t(:);

f=10e6;

f = Df_fft*round(f/Df_fft);

y = cos(2*pi*f*t);  % I think i need some scaling factor here for  +1.5 V --> 0 dBFS

if n_en==1
    y_noise = awgn(y,SNR,'measured');
else
    y_noise= y;
end

S = 2*fft(y_noise,N)/N;
S = 20*log10(abs(S));
f = (fs/1e6)*linspace(0,1,N);
plot(f,S)

So my questions are what NFFT points to take.?
as that would also determine the noise floor.. ..
what exactly would the noise floor be. ?. how exactly i can see the -120 dBFS/Hz in the graph or what does it means in terms of the graph...

is this the correct way to add noise ?

Thanks

All help would be appreciated.

Please see this post that details most of the noise considerations for an ADC

What are advantages of having higher sampling rate of a signal?

In particular the relationship for the noise floor given a full scale sine wave:

$$SNR = 6.02 \text{ dB/bit} + 1.76 \text{ dB}$$

Using your numbers:

1.5Vpp sinewave. If we assume this sinewave is full scale (any higher voltage would start to add clipping noise), then this is the level of 0 dBFS. In particular the variance of this or the standard deviation is used in the dB ratio computation according to:

$$SNR = 20Log_{10}(\sigma_s/\sigma_n)$$

or

$$SNR = 10Log_{10}(\sigma_s^2/\sigma_n^2)$$

Since the noise itself is given in units of dBFS, there is nothing to do further with the 1.5V level. This level is important to determine if the signal is indeed at full scale, and if we wanted to determine the noise level in units of volts.

120 dBFs/Hz with a 200 MHz Sampling Rate:

The quantization noise can be reasonably assumed to be a uniformly distributed white noise (see link above). So the total noise in 200 MHz of bandwidth would be:

$$-120 \text{ dBFS} + 10Log_{10}(200e6) = -37.0 \text{ dB}$$

(Note that this is the double-sided noise in the spectrum extending from $-f_s/2$ to $+f_s/2$)

From the formula above we see that this would be 6.16 bits of precision. A 6 bit converter with the sinewave right at full scale would be very close the the added noise that was given in the example.

To "add" such noise to the signal, scale the signal to fit peak-peak with a 6 bit system and round to closest integer to create the resulting quantization noise.

So the +/-1.5V peak to peak sine wave would then be scaled to $0$ to $2^6$, or $-2^5$ to $+2^5-1$ if you want to represent a signed signal.  The noise will occur from either rounding or truncating to integer values. (There is no noise difference in rounding or truncation; truncation simply introduces an offset).

As far as the FFT, each bin in the FFT if not windowed has a bandwidth 1 bin wide. So the FFT noise floor will go down according to the number of bins. This is explained further in this post:

How large FFTs can pull signals out of the noise floor?


why we do multiple aliasing of high order harmonics the general rule for aliasing.
If my sampling frequency is Fs=800Mhz
signal frequency =120Mhz.
for the 4th harmonics 480MHz
its 80MHz above nyquist frequency(400Mhz) thus its mirrored 80MhZ back to 320MHz.
What about 1120MHz the closest nyquist frequnecy to it is 2*400 so its supposed to be mirrored
800-320=480 why we do double mirroing for it?
why it also mirroed to 320Mhz?

Aliasing is more like circular wrapping around rather than just mirroring.  Frequencies wrap circularly between -Fs/2 to Fs/2, an 800 Hz range from -400 to 400 at a sample rate of 800 sps.
120 is inside -Fs/2 to Fs/2, so doesn't need to wrap around at all.
But 480 wraps 80 past 400, and shows up at -320, which looks like 320 if you don't care about phase.
And 1120 wraps all the way around the circles range of 800 Hz, plus another 320, so it shows up at 320.
12000 Hz (12 kHz) would wrap around 15 times and show up identical to DC or 0 Hz.


KNN method for different values of K on a handwritten digit data set We're supposed to use KNN (K Nearest Neighbor) method for different values of K to classify a handwritten digit data set.

The problem is, this is the first time I want to do a project like that in Matlab. Could someone please explain the normal steps I should do for such problem? What values of K should I use? I read some articles on KNN method and have a rough idea about how the method works but I need help to do this and to be able to go through the following problem we should do after this.

The k-NN algorithm all by itself is very simple, and can be used on vectors of any dimensionality to classify them. Let me write it down here for the sake of answer completeness:

store all your vectors from the training set
store the class of each vector together with them
when a new, unclassified, vector arrives, compare it to all the vectors from the training set
select the $k$ neighbours that have the smallest distance to the new neighbour
the class assigned to the new vector is the most common class among those $k$ closest vectors

Now, that said, the value of $k$ you choose depends greatly on the size of your training set. E.g. if you have a training set of size $100$, and set $k=100$, the result will always be the same, equal to the most frequent class in your training set. If you set $k=1$, the new vector will be classified as the geometrically closest vector (Euclidean distance) every time.

If you're just supposed to present the results for different $k$ for the sake of discussion and learning, I suggest you start from $k=1$ and test for all the $k$'s sequentially: it will make a nice graph and a nice analysis. Maybe you can skip the even values of $k$ (for $k=2$, if you get two different classes as the closest, you're unable to make a decision).

What's really important is how are you going to shape your descriptor vectors. Each component in your descriptor vector is a feature you describe numerically which should be able to help you make a decision. You also have the curse of dimensionality to worry about: more components your vector has, slower the decision making process is.

You have two approaches here:

Decide carefully on features that are confirmed to be good at the task (do some research), and limit yourself to a certain number.
Work with features you think would be good, but don't worry too much about how good they are. Then, when you have a big number of features, use a method like Principal Component Analysis (PCA) or Fisher's linear discriminant (FLD) to reduce the dimensionality of your descriptor vectors before the k-NN

Lastly, let me give you some ideas about features you might want to use to form your feature vectors. I'm not guaranteeing for any of them, I'm just writing ideas off the top of my head:

Down-sample (shrink) your sample images to $3 \times 3$ or $3 \times 4$ pixel black and white images, and take the 0/1 values of all the 9 or 12 pixels as features. You might want to do a morphological thining or get the skeleton of the pen stroke before you shrink your images.
Instead of taking each pixel from a down-sampled image, you could take only the count of background/foreground pixels in each row and column in the downsampled image, making for a total of 6 or 7 non-binary features as opposed to 9 or 12 binary ones.
Count the number of holes in your (binarized) image. By "hole", I mean a connected component of background pixels not connected to the image border.
Take a simple measure such as describing the elongation of the object, e.g. the horizontal distance between two horizontally furthest pixels divided by the vertical distance between two vertically furthest pixels.
Use the Hu set of invariant image moments (wikipedia, article) It's very old work, but I've used them while making a letter-recognition project during my studies, and for simple uses such as this they can still be helpfull.


Multiplier in digital signal processing? whenever we start working on Discrete Cosine Transform(DCT), Loeffler algorithm a key role. Furthur if we approach hardware implementation of this algorithm, we find methods based on Distributed Arithmetic and Shift-and-Add. Both of these method stress on multiplierless circuits. In both cases multiplication is either performed by shift and adding or using a ROM which stores the look-up tables.

The point I am asking here is why everyone is neglecting multipliers replacing them by adders,shift registers in circuits. We can make a N-bit multiplier using adders and registers. Can someone highlight the negative impact of using multipliers?

Once upon a time, chips (processor ICs) did not have anywhere close to millions (now billions) of transistors on them.

A typical fast multiplier, say made out of a tree of layers of carry-save-adder cells, requires a lot more logic gates (thus transistors) than a simple adder, thus also requiring a lot more power, and with less of them fitting on a particular silicon die area.

Repeated conditional-add-and-shift is slower (in the same transistor technology) but requires a lot less logic gates (thus less transistors) when done sequentially.

If an (ancient) FPGA has tons of small memories (but no multipliers), the memories might be useful as look-up tables that are faster than the equivalent logic equations implemented in fabric connected LUT cells.


How to create an AR filter in Matlab My goal is to replicate the procedure described on pages 15-16 (1461-1462) in this paper, prior to adaptive mixture ICA (AMICA): Overlearning in Marginal Distribution-Based ICA: Analysis and Solutions

Using MATLAB,

How do I estimate the AR-coefficient for a one-tap AR process?
How do I remove the AR-process from the data? (The sample rate is 1000 Hz)
Is there any basic literature, manual, or tutorial that you could recommend?

Answering in reverse order,

For a tutorial, I recommend: (the all pole case is an AR model)

J. Makhoul, "Linear prediction: A tutorial review," in Proceedings of the IEEE, vol. 63, no. 4, pp. 561-580, April 1975.

doi: 10.1109/PROC.1975.9792

Abstract: This paper gives an exposition of linear prediction in the analysis of discrete signals. The signal is modeled as a linear combination of its past values and present and past values of a hypothetical input to a system whose output is the given signal. In the frequency domain, this is equivalent to modeling the signal spectrum by a pole-zero spectrum. The major part of the paper is devoted to all-pole models. The model parameters are obtained by a least squares analysis in the time domain. Two methods result, depending on whether the signal is assumed to be stationary or nonstationary. The same results are then derived in the frequency domain. The resulting spectral matching formulation allows for the modeling of selected portions of a spectrum, for arbitrary spectral shaping in the frequency domain, and for the modeling of continuous as well as discrete spectra. This also leads to a discussion of the advantages and disadvantages of the least squares error criterion. A spectral interpretation is given to the normalized minimum prediction error. Applications of the normalized error are given, including the determination of an "optimal" number of poles. The use of linear prediction in data compression is reviewed. For purposes of transmission, particular attention is given to the quantization and encoding of the reflection (or partial correlation) coefficients. Finally, a brief introduction to pole-zero modeling is given.
keywords: {Control theory;Econometrics;Frequency domain analysis;Information analysis;Predictive models;Rhythm;Signal analysis;Statistics;Time series analysis;Tutorial},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1451722&isnumber=31191

How does one get rid of an AR process?

The inverse filter to an all pole filter, is a FIR filter using the AR coefficients.  Inverse filtering is my guess of what that means.

Question one,
There are a number of AR estimators in Matlab's Signal Processing tool box. You can look at those, or read the above paper and roll your own.  Essentially, you have to estimate the covariance of adjacent pairs of samples, and then use one of the methods detailed in the paper.  The term "one tap" might be confusing, [1 a0] might correspond to one tap or two, but [1] doesn't make sense while [1 a0] does.


How to Concatenate white noise signal and applying bandpass filter matlab I created a white noise signal. I am trying to use the array I have created with the white noise to play in 2 channels. When the left channel I was trying to have 5000 samples and the right will have 0 and when in the right channel it will have 5000 samples and the left will be 0(which means no sound from this channel). I keep getting an error when trying to concatenate the sounds and I think I messed up with the loop because every time it loops it will do 5000 samples and I only need 5000 total with a 10000 sample pause in between the 2 sounds. I also had trouble writing the band-pass filter into the code but it worked when using the FDAToool. Any help would be appreciated.

fs = 44.1; %Sample Frequency khz
L=5000; %Sample length for the random signal
P=10000; %Sample Pause Gap
mu=0;
sigma=2;

Left_signal = [sigma*randn(L,1)+mu 0];
Right_signal = [0 sigma*randn(L,1)+mu];

%Loop thru all samples
for num = 2:L,
%Time
%t = num *(1/Fs);
Left_signal = [Left_signal; sigma*randn(L,1)+mu 0];
end

Pause = zeros(P,2);

for num2 = 2:L,
%time
%t = num2 * (1/Fs);
Right_signal = [Right_signal; 0 sigma*randn(L,1)+mu];
end

Long_signal = [Left_signal; Pause; Right_signal];

If I understand the question correctly this should solve your problem:

L = 5000;
P = 10000;
sigma = 2;
mu = 0;

left = sigma*randn(L,1) + mu;
right = sigma*randn(L,1) + mu;

y = [left zeros(L,1); zeros(P,2); zeros(L,1) right];    % stereo sound matrix


Creating a 50hz comb filter in MATLAB I have a signal in a 1d vector X, and I need to remove 50Hz and its resonant frequencies. I've looked into using fdesign.comb but I'm not sure how to apply that filter on my data (after building it using design).

Alternatively, I think I can also do this using iircomb but again, I'm not sure how to exactly. If someone could step me through this (and how to then apply it to a signal X) it would be  appreciated.

In order to use fdesign.comb you have to do something like the following to specify the parameters of the filter

D = fdesign.comb('notch','N,BW',8,.01);

but then you must generate the filter as you mentioned in your question.

H = design(D); % H is a DFILT

you can then view the filter like this

fvtool(H)

One way to apply the filter to your data is to use the filter function

y=filter(H.Numerator,H.Denominator,inArray)

My favorite way to design a filter I've never used is with fdatool because you can play around with all of your parameters. Here is how I would do it:

Type fdatool at the command prompt and a window pops up
Specify Notching in the dropdown that initially says Differentiator
Inder IIR choose comb
Input your parameters
click Design Filter

Once you play around and it looks as you wish, you can either export the coefficients or generate the code that would create this filter.

To export coefficients

go to file and click export
make note of the variable names and click export
now you can use the filter command with y = filter(num,den,data);

To generate code

go to file and click generate MATLAB code
save the file

You see in this file something like the following

Fs = 48000;  % Sampling Frequency

N     = 10;    % Order
BW    = 1200;  % Bandwidth
Apass = 1;     % Bandwidth Attenuation

[b, a] = iircomb(N, BW/(Fs/2), Apass);
Hd     = dfilt.df2(b, a);

the line with iircomb is what generates the filter coefficients, you can now use y=filter(b,a,data);

Hope this helps.


How to find quantization noise in my simulation? For college we have to make simple exercise in simulink (matlab). (since my professor is not familiar with simulink, he can't help).

The question itself is, how can I simulate the SNR when uniform quantization is used. To illustrate this:

Because of this quantization there is noise. I need to display this noise. How do i do this?

Additional information: 10 bit quantization. The source has to be 40db.

I know i can get the quantization error when I subtract the output from the input but i do not know how to get the quantization noise.

The quantization error is the quantization noise. It is two different terms for the same thing.


Fail to separate sound signals by FastICA on real-world recording I have written a program to perform FastICA on a stereo WAV file using the code on Python MDP FastICA Example

With the audio examples I get very good results.

Then I try to do real world recording using two computer mono microphones connected to the stereo mic in of my pc by connecting mic 1 to L channel and mic 2 to R channel. I test by playing some music at the background while I am talking in a quiet room.

However, running FastICA does not separate the signals at all. Is it possible that the quality of microphones is too poor? Do I need to do anything to the recorded WAV file (16 bits, signed PCM, 44100Hz) before running FastICA?

You can download the recording here.

ICA in raw form is only suitable for use with phase synchronised observation mixtures. Using microphones as you have described will introduce a phase delay as pointed out by other posters. However this phase delay can be used to great avail. The best known algorithm that deals with stereo separation in the presence of delays is DUET. The links are broken but the references you are looking for are here >http://eleceng.ucd.ie/~srickard/bss.html.

This is the paper you should look for >
A. Jourjine, S. Rickard, and O. Yilmaz, Blind Separation of Disjoint Orthogonal Signals: Demixing N Sources from 2 Mixtures, IEEE Conference on Acoustics, Speech, and Signal Processing (ICASSP2000), Volume 5, Pages 2985-2988, Istanbul, Turkey, June 2000


SC-FDMA in LTE system I am trying to understand how SC-FDMA is used in LTE system. Every description I found so far refers to DFT and IDFT.

However, when I checked 3GPP TS 36.211 chapter 5.4 (SC-FDMA baseband signal generation) and 6.12 (OFDM baseband signal generation) I found that there is hardly any difference between given equations (I did not copy the equations here because they are quite long) - OFDMA takes into consideration DC subcarrier, in case of SC-FDMA term "k + 1/2" is used wheres OFDMA uses simply term "k" (k can be denoted as number of subcarrier). Moreover the document does not say anything about neither DFT nor IDFT.

In case of OFDMA, I can imagine how the given equation can be considered as IDFT. But what about SC-FDMA? It looks like that modulated symbol is mapped in between two subcarrier. Can somebody explain me this issue in details?

Edit:

I was interested in how DFT and IDFT can be mapped onto 3GPP TS 36.211. In other words, which part of the standard implies usages of DFT and IDFT. Based on the answer given by AlexTP I understand that DFT usage is implied by chapter 5.3.3, whereas IDFT is described in chapter 6.12 (I am refering to version 13.4 of the standard). My first guess was that both DFT and IDFT are implied by chapter 6.12. Extra information that only PUSCH uses SC-FDMA was very useful.

There are some subtleties though. However, it is out of the scope of this question.

First of all, in refering to standard specification, please mention the version because the numeration between versions can be different.

My answer should be only a comment but I think with images it would be easier for you. Everything in my answer is for 36.211 v8.7.0.

The section 5.6 SC-FDMA baseband signal generation and 6.12 OFDM baseband signal generation are IDFT. The difference is the values of $\left\lbrace a_{k,n} \right\rbrace$. Ignoring MIMO stuffs in DL, in OFDM case, $\left\lbrace a_{k,n} \right\rbrace$ are mapped directly to subcarrier. In SC-FDMA case, $\left\lbrace a_{k,n} \right\rbrace$ are outputs of a supplementary block Transform precoder as in the image below.

This transform precoder is implemented as a DFT and its formula can be found in the section 5.3.3 of 36.211 v8.7.0. It is specified only for PUSCH data (if my memory is correct, among UL physical channels, only PUSCH data uses frequency-generated SC-FDMA).

I hope my answer covers what have been asked in your question. For more details (k+1/2, null DC, ...) I suggest a very good book "LTE – The UMTS Long Term Evolution From Theory to Practice" which is surprisingly can be found by a simple Google search.


Butterworth lowpass filter In order to understand the behaviour of the butterworth (lowpass) filter I generated a "spectrum" consisting of 1's. I then went to time domain (ifft), applied the filter and transformed back (fft). I chose the window (maximum) frequency to be one. So i expected that the windowed spectrum should be equal to the actual spectrum of ones.
Instead I got almost the spectrum, but at some point I got a huge line (I expect it to be spectral leakage?!). Also, the values are slightly moved up... Unfortunately I haven't understood how to solve this problem.

Here a sample source code from which I started:

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt

wavelength = range(10000)
spectrum = [1]*len(wavelength)
filter_order = 3 #2
Wn = 1
inter = np.fft.ifft(spectrum)
b,a = signal.butter(filter_order, Wn, btype='low', analog=False)
inter_w = signal.filtfilt(b, a, inter, method='gust') # filtered
spec_w = np.fft.fft(inter_w) # filtered spectrum

plt.plot(wavelength, spectrum, label='spec')
plt.plot(wavelength, spec_w, label='spec_w')
plt.legend()
plt.show()

And here a small plot. As stated I don't understand, why "spec_w" is (i) bigger than one and (ii) there is a peak that shouldn't be there.

Filter_order 3:

Filter_order 2:

Everything you see are numerical artefacts. A "low pass" filter with Nyquist as cut-off frequency is simply a connection. However, you're asking for a recursive filter of order $2$ or $3$, so what you get is a filter with identical numerator and denominator polynomials (i.e., pole/zero cancellation). Check your filter coefficients, a and b should be identical.

So what you have is an unnecessarily complicated implementation of the identity operation, which results in numerical errors. The numerical errors are caused by the "filtering" operation (made worse by the double filtering implemented in filtfilt), and by the combination of FFT and IFFT.

Also note that you should plot the absolute value of spec_w, because due to numerical inaccuracies spec_w might have a small imaginary part and I don't know what your plot function does with a complex-valued vector.


Normalization purpose in signal processing What is the purpose of normalizing the signal?

If we have two signals on hand, how is it used when comparing these two signals?

Normalization is basically bringing the two signals to the same range or a predefined range. A typical example of a predefined range is the statistical perception of the normalization, which is transforming the signal so that its mean is $0$ and standard deviation is $1$. After such transform the canonical form is obtained. In layman's terms, such the operation consists of subtracting the minimum value and dividing by the range (difference between the maximum and minimum values).

Transforming all signals to such canonical form eases and robustifies the process of comparisons as well as serving too different needs such as visualization, and analysis. If we are to speak in terms of image normalization the effect would be the following:

Un-normalized:

Normalized:

Note that, the distinction between the colors is clearer and more apparent.


resampling from fractional to integer sampling rate I have data recorded at 1000Hz and other physio data recorded, somehow, at 1960.8 Hz (510 microsecond samples). Ideally, I would like to be able to downsample the 1960.8 hz signal to line up with the 1000Hz signal.  Any suggestions? I am familiar with matlab or Python.

You can use the resample command mentioned as follows to accomplish this:

y=resample(x,1250,2541); % x is the input signal


How to recover signal repeated throughout other signal? I have an audio signal that contains, at various time offsets, various other short audio signals, often repeated, in addition to other audio content. Here is an example of how that could look like:

   0  1  2  3  4  5  6  7  8  9
A: x        x              x
B: x                       x
C:       x  x

This terrible representation shows that A and B appear at 0 seconds, C appears at 2 seconds, A and C appear at 3 seconds and so on. These signals will additionally be scaled at every occurrence (their amplitude could be different at each occurrence) and possibly have an amplitude envelope applied to them.

Is there an algorithm to recover that takes all the position of the occurrences and recovers (to the best of its ability) A, B, and C?

To be more explicit in the task I have to solve, I have an old electronic music track I produced, but I lost the project file. Like most electronic music, it's a mix of synthesized (pseudoperiodic) sounds, electronic drums hits (called, confusingly, samples), vocals, effects and maybe some acoustic instrument. All those signals get summed together and processed in a "mastering chain": the thing relevant to this problem is that dynamic processing is applied to increase the average loudness, and saturation (in the case of this track, just soft clipping) is added to squeeze another bit of volume.

I need to change one of the drum hits. I don't need to recover it perfectly, because even if a bit of the old drum hit and a bit of noise is still there, the new drum hit that I'm going to replace it with will cover it.

This is a source separation problem, which, in general is quite difficult to solve. One way to attack it is to use non-negative matrix factorization (NMF) with short-time Fourier analysis. Please see this question and a sample script here.

The overarching idea is to try to isolate your component sounds A, B, C,... using a time-frequency decomposition (or, perhaps, some other kind of subspace analysis using a different transform, like wavelet transform).

If the components don't always all play at exactly the same time, or have different onsets, you can hope to separate them temporally. If they don't all have identical frequency content, you have some hope of separating them in frequency. Or, a combination of both, which is what NMF of the STFT achieves.


Calculating normalized edge frequencies for filter design in Matlab What is the proper formula for converting ordinary frequency to normalized edge frequency in matlab for use in filter design commands like cheby2 etc
The wikipedia page shows three different formulas for normalized frequency,i higlighted all those in attached snapshot
And Matlab magnitude plot has written  (x pi rad/sample) below plot,and i highlighted that also.
So which one of three formulas of wikipedia ,should we use in Matlab?

This is from the mathworks documentation on cheby2

For digital filters, the stopband edge frequencies must lie between 0 and 1, where 1 corresponds to the Nyquist rate—half the sample rate or $\pi$ rad/sample.

For discrete-time signals, we use the normalized frequency $\omega$ in radians (per sample), defined as
$$\omega=\frac{2\pi f}{f_s}\tag{1}$$
where $f$ is the actual frequency in Hertz, and $f_s$ is the sampling frequency.
Matlab uses $(1)$ normalized by $\pi$, i.e., edge frequencies etc. are defined by W$=2f/f_s$. The value W$=1$ corresponds to $f=f_s/2$, which is the Nyquist frequency.
Example: The cut-off frequency is $f_c=500$ Hz and the sampling frequency is $f_s=2000$ Hz (samples/second). According to $(1)$, the corresponding normalized frequency in radians (per sample) is
$$\omega_c=\frac{2\pi f_c}{f_s}=\frac{\pi}{2}\tag{2}$$
Consequently, when calling a Matlab routine you would use Wc$=\omega_c/\pi=\frac12$.