texts:
How should 4 point polyBLAMP work. It does not really seem to dampen aliasing.

I have an floating point digital audio oscillator that makes a +1.0f jump at the end of every period.

My question is, have I misunderstood how it should work?

I have been using these as reference:
http://www.dafx17.eca.ed.ac.uk/papers/DAFx17_paper_100.pdf
https://www.researchgate.net/publication/307990687_Rounding_Corners_with_BLAMP

Here is how I thought it should be done: (the code gets run every sample)

float t = phase / period; //normalized position in period
float dt = deltaPhase / period; // normalized phasestep

float polyBLAMP = 0.0f;
if (t < dt) {
    float d = (t)/dt; // fractional delay
    polyBLAMP = (3.0f*powf(d,5.0f)-10.0f*powf(d,4.0f)+40.0f*powf(d,2.0f)-60.0f*d/2.0f+28.0f)/120.0f;
} else if (t < dt*2.0f) {
    float d = (t-dt)/dt;
    polyBLAMP = (-powf(d,5.0f)+5.0f*powf(d,4.0f)-10.0f*powf(d,3.0f)+10.0f*powf(d,2.0f)-5.0f*d+1.0f)/120.0f;
} else if (t > 1.0f - dt) {
    float d = (dt-1.0f+t)/dt;
    polyBLAMP = (-3.0f*powf(d,5.0f)+5.0f*powf(d,4.0f)+10.0f*powf(d,3.0f)+10.0f*powf(d,2.0f)+5.0f*d+1.0f)/120.0f;
} else if (t > 1.0f - dt*2.0f) {
    float d = (2.0f*dt-1.0f+t)/dt;
    polyBLAMP = powf(d,5.0f)/120.0f;
}
float output = osc(phase)+polyBLAMP;


texts: Does Forward Error Correction (FEC) affect upsampling factor?

texts:
Is the redundant data introduced during forward error correction (FEC) counted as upsampling?

I'm taking a binary data stream, doing convolutional encoding of rate 1/2 on it, then making it NRZ signal and then upsampling it. In this page, the relation between $EsN0$ and SNR uses a two variables $T_{sym}$ and $T_{samp}$ as $T_{sym}~/~T_{samp}$ which indicates the upsampling factor. So should I include the rate 1/2 of FEC in $T_{sym}~/~T_{samp}$  calculation or not?


texts: Why does stability only consider, whether the system is bounded? That doesn't mean the response is smooth?

texts:
Why does (BIBO) stability only consider, whether the system is bounded? That doesn't mean the response is smooth?

That is, one could have system whose values are bounded. Still it could exhibit "chaotic" behaviour through its parameters.

So what is this kind of stability/instability concept really helpful for?


texts: I have a 3D shape in a 3D binary image. Therefore, I have a list of all of the x,y,z points. If I am to analyze this shape for various identification, such as "sphericity", volume, surface area, etc., what are some of the choices do I have here?


texts: Using Kalman Filter to fuse two sensor readings of the same type

texts:
I have always used Kalman Filter to smooth a signal comes from one sensor only.

I would like to know if Kalman Filter could be used to fuse data coming from two different sensors that provide the same type of reading. An example is fusing the position data that comes from the GPS with the position data that is calculated from the IMU. In the light of this example, I wish to have the following:

State vector: [x;y]
measurements vector: [x_imu;y_imu;x_gps;y_gps]

Let me know if you need more information/clarification.


texts: Most accurate way to find RMS of signals fundamental

texts:
I would like to know which algorithmn is the most accurate one in order to calculate the fundamental RMS of two signals (f. e. voltage and current). Is it the standard FFT approach with a rectangular window or may there be some more precise algorithmn?

In this example, the resulting RMS signals may be used further for calculating the signal of the apparent power. The following Python code may help to illustrate the question:

import numpy as np

#sampling frequency in Hz
Fs = 16000

#fundamental frequency in Hz
F1 = 50

#number of samples
samples = 1000000

#generate voltage and current signals
t = np.arange(samples / Fs)
u = np.sin(F1 * 2.0 * np.pi * t) + 0.5 * np.sin(80.0 * 2.0 * np.pi * t) + 0.1 * np.sin(10.0 * 2.0 * np.pi * t) #signal may have additional harmonic and sub-harmonic content
i = np.sin(F1 * 2.0 * np.pi * t) + 0.2 * np.sin(120.0 * 2.0 * np.pi * t) + 0.3 * np.sin(30.0 * 2.0 * np.pi * t) #signal may have additional harmonic and sub-harmonic content

#apply some algorithm to find accurate RMS of the signals fundamental
#t_rms_50Hz = ...
#U_rms_50Hz = ...
#I_rms_50Hz = ...

#get apparent power S
S_50Hz = U_rms_50Hz * I_rms_50Hz


texts: Shape of Kaiser windowing function FFT response

texts: I'm trying to determine the shape of the FFT response of a signal which had a kaiser windowing function applied to it. What I mean by this is if the shape is a normal bell curve or a Lorentzian distribution shape. I want to be able do a curve fit to the FFT response of the windowed signal for an analysis, but cannot find information on the shape of it. Here is an example of what the FFT response looks like:


texts:
I saw some questions (and answers) on this subject, but they were all about a specific example and I'm not sure I understood.

I'm trying to understand the meaning of computing the DFT of an array representing a filter (let's say 1D).

I thought I need to think of the array (let's say [1,2,1]) as the samples of the function (signal), and just compute the DFT. But I'm not sure any more.

let's say I have:

cos(2pi(4x+4y)\256)

so it's DFT is some a at (4,4) and (-4,-4).
now I want to compute g*[1,2,1].
so I want to compute the DFT of f = [1,2,1], (by the convolution theorem).
I used the formula, but couldn't see the connection to G, the index +-(4,4) etc.

also, if f=[1,...,1] (all ones), I thought I should get 1 at F(0), and zero anywhere else. as it is a fixed function.
but this is not the case. can you explain why?

thank you!


texts: Relationship between interpolation accuracy, impulse response and frequency response

texts:
Global interpolation or sinc interpolation is an ideal filter since its frequency response is a rect function.
The impulse response of this filter is the sinc function (same as the coefficients of the interpolator).

Local interpolation methods are constructed by applying some window to the sinc function.
For example, according to https://ccrma.stanford.edu/~jos/Interpolation/Relation_Lagrange_Interpolation_Windowed.html, for uniformly spaced samples and finite $N$ , Lagrange interpolaton is equivalent to windowed sinc interpolation using a binomial window.
The impulse response of this filter are the coefficients of the interpolator.

I am interested to know what is the relationship between the impulse response, accuracy of the interpolation and frequency response of the filter/local interpolator, especifically for Lagrange interpolation.
I cannot see this clearly.


texts: Impulse Response of real coefficient, LTI System

texts:
I'm trying to obtain the impulse response $h[n]$ of a system whose frequency response is $H(e^{j\omega})=R(\omega)e^{-25j\omega}$. I believed that $h[n]=h[n-25]$, would be the correct answer, however I was told it would be something along the lines of $h[n]= h[50-n], 0<=n<=25$.
Can somebody explain why?


texts:
I am solving the following problem as an exercise.
I have a signal
$$z(t) = x(t)\cos(wt) + y(t)\cos(ut)$$
and we assume $x,y$ to be band-limited continuous time signals. In the frequency domain we have that both signals are $0$ if the frequency $\omega \geq \omega_{max}$. We try to reconstruct $x$ by passing $r$ through a bandpass filter $B(\omega)$ that is $1$ if $\omega$ is $\omega = w$ or $\omega = -w$ (it actually is $1$ in an interval centered around $w$ and $-w$ but I didn't know how to write it since I have no information about the width of the interval. The outcome is then multiplied by $\cos(wt)$ and then passed through a lowpass filter $L(\omega)$ which is $2$ if $-\omega_{max}\leq \omega\leq  \omega_{max}$. What do we need to assume on $w,u$ in order for $x$ to be reconstructed?
In the frequency domain, I have that
$$Z(\omega)= \frac12\big(X(\omega+w)+X(\omega-w)+Y(\omega-u)+Y(\omega+u)\big)$$
So picotrially I have 2 copies of $X$ centred at $w$ and $-w$ and two copies of $Y$ centred at $-u$ and $u$ with half the amplitude, right?
Now, the first bandpass filter should give me only the two copies of $X$, since it is centred around $w$, then there is this multiplication with $\cos(wt)$ again. What does this do? Doesn't it simply make yet another two copies of $X$ centred around $w$? Does it give me back the original amplitude getting rid of the $1/2$? Then, how can a lowpass filter allow me to reconstruct the signal when I have $2$ copies of it?


texts: How do I find a 50Hz sin/cos signal in a measurement among other frequencies and noise?

texts:
I measure magnetic fields with high precision, in a gradiometer setup with fluxgate magnetometers, sampling at 2000Hz. I can have big data sets if I want to - one is 600000+ samples long. I need to find the phase and amplitude of a 50Hz sine signal. However there are other (low frequency - 2.145Hz) stronger signals superimposed, along with noise. Because of those low frequency signals the DC offset is hard to eliminate, too. As a result, an FFT gives no results yet.

Now it occured to me that it might be easier to find the 50Hz signal I look for (which is several orders of magnitude smaller then the DC offset and the other 2.1Hz signal) and determine its phase and amplitude, in stead of filtering out the rest.

What are good approaches and algorithms for that? I code this in python.

edit: here is my python (with numpy) code for my carlos loop algorithm.

   def decode_with_carlos_loop(self, frequency, periods_num):
        offsetfree_data = self.file_data
        shape = offsetfree_data.shape
        sampling_freq = 2000

        chunk_size = int(sampling_freq/frequency)
        alpha_step = .075
        indices = np.arange(chunk_size)

        for sensor_cnt in range(0, 4):
            for axis_cnt in range(1,3):
                error_sin = error_cos = 0
                cos_lo = sin_lo = []
                theta_last = 0
                theta = np.zeros(int(shape[1] / chunk_size) + 1)

                for chunk_cnt, chunk_pointer in enumerate(range(0,shape[1],chunk_size)):
                    chunk = offsetfree_data[sensor_cnt, chunk_pointer:chunk_pointer+chunk_size, axis_cnt]

                    cos_lo = np.cos(2*np.pi*frequency*indices/sampling_freq + theta[chunk_cnt])
                    sin_lo = np.sin(2*np.pi*frequency*indices/sampling_freq + theta[chunk_cnt])

                    # mean is the lowpass filter
                    error_sin = np.mean(np.multiply(chunk, sin_lo))
                    error_cos = np.mean(np.multiply(chunk, cos_lo))
                    error = error_cos*error_sin
                    theta_last = theta[chunk_cnt+1] = theta[chunk_cnt] - alpha_step*error
                print("sensor:axis:" + str(sensor_cnt) + ":" + str(axis_cnt) + "\tlast theta: " + str((theta_last + np.pi) % (2 * np.pi) - np.pi) + "\tcos_lo: " + str(cos_lo[0]) + "\tsin_lo: " + str(sin_lo[0]) )


texts: Why are some LDPC codes more immune to noise than others?

texts:
When one analyzes an LDPC code, it is usual to plot the BER versus the $E_b / N_0$. My questions are:

Why are some LDPC codes better than others in terms of noise immunity? After all, they are just defined by matrices $\mathbf H$.
What properties or parameters of this parity-check matrix have influence on the BER for a given SNR?


texts: Is this two input discrete system linear?

texts:
Given two inputs $\: x_1[n]\: x_2[n]\:$Is the system  $\:y[n]=x_1[n]\times x_2[n]\:$ linear ?
My Approach:
$(x_1\times x_2)[n]=S_1[n]\rightarrow Y_1[n]$
$(x_3\times x_4)[n]=S_2[n]\rightarrow Y_2[n]$
$(ax_1\times x_2)[n]+(b x_3\times x_4)[n]= a  S_1[n]+ b S_2 [n]  \rightarrow a Y_1[n]+bY_2[n]$
If that is correct it is linear. But I am not sure about this solution. I will appreciate any help.


texts: Estimate input from measuring only RC filtered input

texts:
I am interested in knowing the step amplitude on a signal which can only be measured through a hardware RC filter. The measured signal follows this formula:

$$
V(t) = V_0 \left(1 - e^{-\frac{t}{RC}}\right)
$$

Is it possible to estimate the step amplitude, $V_0$, given that we can perform two measurements on the almost linear part of the rising edge (with a known delay between them)?


texts:
Sorry if this is a stupid question but so far i have learned that after getting the results of a FFT (an array of complex numbers) you can get the frequency of a sound sample via the following:

1) Calculate the peak for each complex number: $\sqrt{\mathrm{real}^2 + \mathrm{imag}^2}$

2) Get the index of the array element with the highest peak.

3) Index $\times$ mhz of the sample audio $/$ number of elements in the array

This seems to work however I'm wondering if the second-highest peak would be the second "most present" frequency, and so on. So for example, if I had a guitar chord, and I got the top 6 frequencies, would those be the notes of a 6 stringed chord?


texts: Calculate sum of absolute difference per pixel

texts:
I am trying to implement a research paper just to verify my results. I have shared the code (written by me) and the description of the process given in the research paper. Unfortunately my results are a bit different from his results that means I am doing something wrong. Kindly guide me in the right direction.

Here is the link to the paper: A parametric model for perceptual video quality estimation.

I have tried for mobile video sequence which can be downloaded from the following link: Yuv video sequences.

He gets 3.600 as SAD/pixel value while I am getting around 1.

Code:

fileName = 'mobile_cif.yuv';
    width=352;
    height=288;
    numFrames=30;
    SAD=0;

    %for loop to traverse & process from frame '1' to 'last' frames
    for t = 1 : numFrames
    if t < numFrames %(Need to traverse to n-1 frames)
    result=loadFileYuv(fileName,width,height,t);
    currFrame = result.cdata;    %reading Current frame
    result=loadFileYuv(fileName,width,height,t+1);
    nextFrame = result.cdata;    %reading Next frame
    % First, take the absolute value of the difference at each pixel
    myAbsDiff = abs(double(currFrame) - double(nextFrame));
    % Then, sum over all pixels
    out = sum(myAbsDiff(:));
    out=out/(width*height); %(Taking Average with respect to pixel)
    out=out/64;%(Dividing By Block size to normalize)
    disp(out);
    SAD=SAD+out;
    else
        continue;
    end

    end
    disp(SAD/numFrames);  %(Taking Average with respect to Frames)

Process Description

Sum of absolute difference (SAD) is a simple video quality metric used for block comparison and for moving vectors calculations. Each frame is divided into small blocks (i.e. $8 \times 8$ pixels) and for every block in one frame the most similar (minimum SAD) block in next frame is find. This minimum sum of absolutes differences is assigned as the SAD for each block in each frame (up to the $n − 1$ frame). Then all the SAD values are averaged for each frame and for all the frames in the clip, and divided by the block area, for normalization.


texts:
I'm learning Kalman Filter for a week now. I just discovered that EKF (extended Kalman Filter) might be more appropriate for my case.

Le't suppose I'm applying KF/EKF for variometer (the device that tells planes and parachuters what's their vertical position and velocity). In my case I've generated some sample data: the first few seconds he (the parachuter for e.g.) is falling (the velocity is positive) then he is rising (velocity is negative).

As far as I can tell this system is linear. So should I use KF or EKF?


texts: Can the j be equal to k in the first formula ?In MISO system

texts:
In this system model,the received signal available for RF EH at each RX_k is given by

However,its SINR formula is given by

We can find that author said the j can't be equal to k,but their range are the same in the second formula .However,in the first formula,i am not sure that is the j equal to the j in the second formula?
I mean,now i know j which is in the second formula can't be the same as k,but how about j in the first formula?can it be equal to k ?

Well,i think j can be equal to k in the first formula,however,i found that if i write the code in "j can be equal to k" this situation,i can't get the optimal value every time,sometimes,the optimal value will be "Inf".however,if i write the code in "j can't be equal to k" this situation,i can get the optimal value every time,it is strange,so i need someone to tell me can j be equal to k in the first formula?
Paper link: https://arxiv.org/pdf/1805.08898.pdf


texts: How do you handle negative frequencies in the power spectrum of a complex signal?

texts:
When we apply the DFT operation on a real signal $x[n]$ to get $X[k]$,
then take the square magnitude of the $X[k]$, $\lvert X[k]\rvert^2$, the power spectrum is symmetrical.
You can take the positive frequencies or negative frequencies as the frequency information in $X[k]$.

However this is not true for complex valued signals; the power spectrum is not symmetrical.

In this case, how would you determine the frequency components in original signal?
Can we just drop the negative frequency part?


texts: Unexpected imaginary part in the fft of a zero-padded cosine

texts:
I simulated a cosine waveform $ y = \cos(\omega t) $ and I applied the FFT algorithm to it. As expected, I have a frequency peak at $\pm \omega$ only in the real part, and nothing in the imaginary part. As expected, if I zero-pad the array, the peak in the real part becomes a $\mathrm{sinc}(\omega D)$ (where D is the length of the window). However, in the imaginary part, a sort of anti-symmetric $\mathrm{sinc}$ appears. I don't understand why, because I was expecting the imaginary part to stay zero since the Fourier transform of the cosine is purely real.
Why is this happening?
Here a plot of the imaginary and real parts of the FFT of the described signal. In orange: FFT of the original signal. In blue: FFT of the zero-padded signal.

Here is my code, written in MATLAB language.
dt = 0.02;
fs = 1/dt;

start = 0;
stop = 1;
t = start:dt:stop-dt;

w = 2*pi*10;
y = cos(w*t);

N = 10000;
y_pad = [y zeros(1,N)];

Y_PAD = fftshift(fft(y_pad))./length(y);
df_pad = fs / length(y_pad);
f_pad = (0:(length(y_pad)-1))*df_pad;
f_pad(f_pad >= fs/2) = f_pad(f_pad >= fs/2) - fs;
f_pad = fftshift(f_pad);

Y = fftshift(fft(y))./length(y);
df = fs / length(y);
f = (0:(length(y)-1))*df;
f(f >= fs/2) = f(f >= fs/2) - fs;
f = fftshift(f);

figure;
subplot(2,1,1);

plot(f_pad, real(Y_PAD));
hold on
plot(f, real(Y));

ylabel('Re(fft)')
xlabel('f (Hz)')

subplot(2,1,2);

plot(f_pad, imag(Y_PAD));
hold on
plot(f, imag(Y));

ylabel('Im(fft)')
xlabel('f (Hz)')

Edit: By zero-padding only to the right of my signal I am breaking the even parity of the cosine and $y[N-n]=y[n]$ does not hold anymore (see Hilmar's answer). Therefore, I tried to zero pad both to the left and to the right, in order to preserve the parity
N = 10000;
y_pad = [zeros(1,N) y zeros(1,N)];

The result is the following and it looks weird: the imaginary part is still not zero and the real part oscillates between positive and negative values. Why this zero-padding does not work as expected?


texts: Does the auto-correlation function of stationary random process always converge?

texts:
The auto-correlation function of the stationary random process only depends on the time difference $\tau$.

http://web.ntpu.edu.tw/~yshan/chapter6_han.pdf

64th slide of this lecture note mentions that for a zero-mean process, the autocorrelation converges to zero as $\tau$ goes to infinity.

I want to know why.


texts:
I've found a material explaining the sample and hold operation pretty well, but there's one thing I can't understand.

http://www.springer.com/cda/content/document/cda_downloaddocument/9781848001183-c1.pdf?SGWID=0-0-45-484012-p173780368

Next, we would like to take this ideally sampled signal and hold the values
between samples. As we have a train of impulses with the correct areas, we
need a “block” that takes an impulse with area A, transforms it into a rectangular
pulse of height A that starts at the time at which the delta function is
input to the block, and persists for exactly Ts seconds. A little bit of thought
shows that what we need is a linear, time-invariant (LTI) filter whose impulse
response, h(t), is 1 between t = 0 and t = Ts and is zero elsewhere.

How do we know this is the impulse response that will give us a pulse of HEIGHT A, not area A? Could anyone explain this "a little bit of thought" part? Thanks in advance.


texts: Finding the maximum of a real valued trigonometric polynomial

texts:
I have been posting all over SE looking for a good method of finding the maximum of a trigonometric polynomial. A trigonometric polynomial is the sum of a set of sinusoids,

$$a_0 + \sum_{k = 1}^{N} a_k \cos(k \theta + \phi_k)$$

and can also be written in the form

$$ \sum_{k=-N}^{N} c_k \, e^{i k \theta} $$
where $c_k$ = $\bar{c}_{-k}$

Substituting $z = e^{i \theta}$ turns this into the typical form of a complex valued polynomial.

The method I am currently using is to take the find the roots of the derivative of the polynomial, then evaluate the original polynomial at these locations to find the maximum. While quite fast with MATLABs roots command, it is not quite fast enough.

I need a method where I can trade accuracy for speed. Does anyone know of any other fast methods? This is a similar problem to estimating the crest factor of a signal. In my problem $N$ ranges from 3 to 15.


texts:
While learning sampling theory - I noticed that examples of continuous signal sampling always achieved the goal via multiplying the signal with a "Dirac Comb".
I was intrigued by the requirement to use a Dirac Comb instead of a "Ones Comb" - I.E: multiply the signal by the value '1' at the sampling times.

The explanation I got: "That's necessary because multiplying by a Dirac function is equivalent to Convolution in the frequency domain which maintains the spectrum of the original sampled signal".

This made a lot of sense to me...But than, when I proceeded to learn about discrete domain sampling - I was surprised to find out that examples used a "Ones Comb" to sample a signal (not the "Dirac Comb").

Why the difference ? How does discrete time sampling get away with multiplying by a train of discrete ones ?


texts: Given Gain function, how to design a causal, stable and minimum phase IIR filter?

texts: I am given $|H(\omega)|$, I wonder if minimum phase stable causal filter is unique and how to calculate it.


texts: Least Squares and Auto and Cross Correlation

texts:
I am trying to understand why auto and cross correlation helps find the best fit line in least squares. I have an equation as stated here:

$Ax=b$ -- I have not exact solution, so I use the least squares equation.
$$A^TAx=A^Tb$$
I can solve and use the equation but I am trying to figure out why cross and auto correlation helps. Can anyone shed some light on that? For example, what do cross and auto correlation represent?

I am learning about signal processing and linerisation of a signal. $A$ is the input unpredistoted signal and $b$ is feedback distorted signal. $x$ is the coefficients that will be used to weight the Look Up Table. $Ax=B$ is solved in the adaptive algorithm box, in the link below.

http://www.altera.com/end-markets/wireless/advanced-dsp/predistortion/wir-digital-predistortion.html


texts:
I have the following image I and I want to add additive noise using the array N below. Can you please give me a calculation example of how to achieve this?

Do I have to move from image space to frequency space?

Thank you in advance!


texts: Why do waveforms that are symmetrical above and below their horizontal centerlines contain no even-numbered harmonics?

texts: All About Circuits site states that waveforms that are symmetrical above and below their horizontal centerlines contain no even-numbered harmonics. Can somebody explain this mathematically, or point to a resource? I do not know much about dsp, but understand what a Fourier transform is.


texts: What is the physical meaning of the convolution of two signals?

texts: If we convolve 2 signals we get a third signal. What does this third signal represent in relation to the input signals?


texts: Bluestein's algorithm to evaluate the DFT from $f_o$ to $f_o + k\Delta_F$

texts:
Briefly, the convolution between $x(nT) e^{-j2\pi f_o nT} e^{-j \pi \Delta_F Tn^2}$ and $c(nT) = e^{j \pi \Delta_F T n^2}$ multiplied $e^{j \pi \Delta_F T k^2}$ allows me to find the DFT $X(f_k = f_o + k\Delta_F)$.
Let $N_x$ be the number of samples of $x(nT)$, I found that $c(nT)$ must be at least of length $2N_x$ but I don't know why. Can someone tell me how do I choose the number of samples of $c(nT)$?


texts: Usefulness of Matrix Notation for Linear Periodically Time Variant Transformations

texts:
I am quit puzzled about a notation I found in Gardner's 1986 book titled Introduction to Random Processes. It is in the chapter on Cyclostationary Processes at section 12.4 which pertains to Linear Periodically Time-Variant Transformations:

  (...) many modulation systems can be modelled as the scalar response
  of a multi-input LPTV transformation with stationary excitation. This
  includes amplitude modulation, phase and frequency modulation,
  quadrature amplitude modulation, pulse-amplitude modulation, all
  synchronous digital modulations such as phase-shift keying,
  frequency-shift keying, and so on. Consequently, the study of
  cyclostationarity is facilitated by general formulas that describes
  cyclic spectra in terms of the parameters of LPTV transformations.
  This includes cyclic spectra that are generated by LPTV
  transformations of processes that exhibit no cyclostationarity, as
  well as cyclic spectra that are transformed by LPTV transformations of
  processes that do exhibit cyclostationarity. Let us consider the LPTV
  transformation $$ Y(t) = \int_{-\infty}^\infty
\textbf{h}(t,u)\textbf{X}(u)du $$ for which $\textbf{X}(t)$ is a
  (column) vector excitation, $Y(t)$ is a scalar response, and
  $\textbf{h}(t, u)=\textbf{h}(t+T, u+T)$ is the (row) vector of
  impulse-response functions that specify the transformation.

The subject of LPTV Transformations is new for me and I don’t understand the usefulness of using non-scalars for $\textbf{h}(t, u)$ and $\textbf{X}(t)$.

Is this a commonly used notation?

Why should a matrix product be used?

I have been trying to make sense of this for some days now and I have really made no progress whatsoever. I hope someone here has an idea...

With thanks!


texts: Separation of a delayed mixtures of signals

texts:
I have three emitters in three different positions that emite three signals $S_1(t), S_2(t), S_3(t)$ with
$$S_1(t)=S_2(t)=S_3(t)$$
Two receivers recive $R_1(t), R_2(t)$ with
\begin{align}
R_1(t)&=r_1(t-D_1)+r_2(t-D_2)+r_3(t-D_3)\\
R_2(t)&=r_1(t-D_1)+r_2(t-D_2)+r_3(t-D_3)\\
r(t)&=\operatorname{conv}\big(S(t),\rm{Imp_r}\big)
\end{align}
Where $\rm Imp_r$ is a realization of the channel impulse response.
My goal is to separate $R(t)$ and get $r_1(t-D_1), r_2(t-D_2), r_3(t-D_3)$ with the exacte delays. ICA and SCA methods works only on instatinous mixtures, in this case we have a delayed mixtures. is there any other methods to solve this problem ?


texts: why 2 vectors are not the same in matlab?

texts:
I have 2 vectors:

a = 1.0e-04*[0.110000000000000   0.140000000000000   0.140000000000000   0.140000000000000];
b =         [0.000011000000000   0.000014000000000   0.000014000000000   0.000014000000000];

by look at by eyes, it seems similar; but when I use this command:

a == b

then I got:

ans =
  1×4 logical array
   0   0   0   0

The expectation is:

ans =
      1×4 logical array
       1   1   1   1

Does anyone knows what wrong with this command in the matlab?

Thank you.


texts: What is a good way to distinguish aircraft noise from other sounds?

texts:
The FAA's NextGen routes, which were modeled, not measured, have created a bunch of aircraft noise for tens of thousands of people. I want to make it easy for those people (including me) to buy a cheap USB microphone, plug it into a laptop or PC, detect aircraft noise, correlate that with publicly available flight-path data, and file complaints with the FAA, airport(s), and legislators. With enough actual (as opposed to theoretical) data, we might be able to get them to roll back the route changes. Or not, but it's worth a try.

To do that, I need some help with the "detect aircraft noise" part of the problem. It's not enough to trigger off noise above a certain dB threshold, because lawnmowers, leaf blowers, trucks, etc. also make noise. I can weed out some over-threshold noise that doesn't last long enough or lasts too long (aircraft take 5-10 seconds to pass overhead, vs. minutes for lawnmowers and a couple of seconds for trucks), but that's a pretty crude approach. And I still need to be able to factor multiple noise sources into individual components.

Correlation would probably be easier with a pair or array of microphones, to establish noise direction that could be tied to flight paths, but that adds a lot more hardware and complexity. Ideally we'd just need one microphone.

Aircraft include airplanes (jets, turboprops, and standard propellor-driven), as well as helicopters.

I'm a good coder (Java expert, former C expert), but have zero experience with DSP or training algorithms. Can you help me get started?


texts: Why is a negative exponent present in Fourier and Laplace transform?

texts:
could anyone explain why there is a need of negative exponent in fourier and laplace transform.I looked through the web but I couldn't get anything.Does anything happen if a positive exponent is placed in these transforms.

While looking through http://1drv.ms/1tbV45S it says that if $s>0$  it becomes a rapidly decreasing function while if $s<0$ it becomes an rapidly increasing functin of t.I couldn't understand that.Can anyone illustrate this.


texts: Can I have some feedback on this research paper?

texts:
I have submitted this paper "Effect of Transducer Positioning in Active Noise Control" to arXiv.

I would like to get some feedback on the methodology of the same. Inputs of any sort are welcome.
Can you please take a look at it and give me some feedback/suggestions/critics?


texts: What happens, when an unlimited bandpass-signal gets mixed?

texts:
In theoretical books there's always a bandpass-limited signal being mixed, that's been low-pass filtered beforehand.

But what happens, if the RF-signal from the antenna isn't band- or lowpass-filtered?

I know for example from decimation (lowpass-filtering and downsampling) that there would be aliasing, which isn't a form of mixing but nonetheless.


texts: I am  a  newcomer in signal processing. I saw that the $L^2$-norm of a signal is also applied as its energy! How is this concept illustrated for those ones who are working in pure math.


texts: Shortest geometric distance from surface in 3d dataset?

texts:
I have a three-dimensional binary image of a collection of discrete, individual voxels ("seeds") contained in a connected 3-dimensional surface ("skin").
(Like a small fruit, with a surface delineated by a one-pixel boundary, that contains seeds.)

The binary matrix was derived from a three-dimensional intensity image of the same size (grayscale).
(It's easy to label the surface voxels "2", the seed voxels "1", and the background "0" using MATLAB's label function).

In MATLAB,
 How do I compute the shortest geometric distance between each individual "seed" and the surface?
In the (crude) 2-dimensional diagram in the link below, each desired distance is represented by the magnitudes of the red lines (in number of voxels).


texts: Detecting changes between two similar images

texts:
I have the following problem:

I need to detect changes between two more or less similar image.

In that case, I want to detect the missing parts.

Lets assume that the displacement is no more that 2px, and there is no significant change in light (but there is some) or camera pose.

What I want as an output, is a binary image about if there is a signifiacant change at that position.

My naive idea was, to apply blur (aroud the size of expected displacement), and compare pixelwisely, but in practice its not working really well. (mainly it has a huge amount of noise around edges)

Any better idea for this problem?


texts: Difference Between Equiripple & Least Squares Design for FIR Digital Filters

texts:
For an efficient and optimized digital FIR filter design, there are two methods available broadly, Equiripple filter design & Least Squares filter design.  A general method for designing  a filter is also Frequency Sampled FIR design but it is not an optimized design

The basic knowledge I have is that Equiripple filter, as the name suggests, has equal ripples in passband & stopband, which means the signal distortion that happens at the edge of the passband due to presence of a large ripple is avoided in Equiripple design BUT, Equiripple design has a large transition band, thus limiting the total passband width.

On the other hand, in a Least Squares design, the transition band width is smaller than for Equiripple design, hence the passband width is more, but the passband ripple are not equi-ripple & exhibit a spike at the passband edge due to Gibbs phenomenon, which causes signal distortion at the edge.

My question is can somebody present or quote the all the differences & advantages over other, in a technical language, for the Equiripple design vs Least squares design of digital low pass FIR filter.


texts: Is autospectrum the same as power spectrum density?

texts:
In my slides about signal processing, there is one that mentions the same thing as the beginning of this and this answer, namely that the Fourier transform of a signal, squared is the power spectrum density of the signal.

In this talk, it is mentioned that coherence is calculated by dividing the squared cross-spectrum by the product of the two auto-spectra.

However, the formula in my slides divides the squared cross-spectrum by the product of the formula of our previously seen PSD and another PSD.

So, is autospectrum the same as PSD? I can find a lot of information about PSD, but not on autospectrum.


texts: How to use values from the output of a block in GNURadio

texts:
I want to use the output values of the top_block as a way of directing my flowdiagram. So I want the flowdiagram to take one path if the output of my processing block is 0 and take a different path if it is non-zero. The if self.blocks_multiply_xx_0.out_sig != 0: statement is the condition I am having trouble with.
I have implemented an if statement within the top block which disconnects and connect parts of the diagram accordingly.

I want the if statement to look at the output of the self.blocks_multiply_xx_0 block and send it to a file sink if it is non-zero, or otherwise send it to a null sink.
I want this to be done for each and every value coming out of the self.blocks_multiply_xx_0 block.

However to do so I need a way of accessing that blocks output value.

How can I go about accessing that value, and assessing which path to take based upon that value, each and every time the self.blocks_multiply_xx_0 block produces an output value?

Is the location of the if statement correct?

Will it check each and every processed sample value or does it have to be moved to a different method?

If so where could it go?
from gnuradio import gr
from gnuradio import audio
from gnuradio import analog
from PyQt4 import Qt
from gnuradio import blocks
import threading
import time
class my_top_block(gr.top_block):
def init(self):
gr.top_block.init(self)
     samp_rate = 32000
     amp = 0.1
     # Blocks
     ##################################################
     self.blocks_multiply_xx_0 = blocks.multiply_vff(1)
     self.audio_sink_0 = audio.sink(samp_rate, '', True)
     self.blocks_file_sink_0 = blocks.file_sink(gr.sizeof_float*1, 'C:\SDR\GNU files\speed_test_SSD\speed_test_SSD.iq', False)
     self.blocks_null_sink_0 = blocks.null_sink(gr.sizeof_float * 1)
     self.analog_sig_source_x_0_0 = analog.sig_source_f(samp_rate, analog.GR_SIN_WAVE, 350, amp)
     self.analog_sig_source_x_0 = analog.sig_source_f(samp_rate, analog.GR_SQR_WAVE, 0.5, 1)

     ##################################################
     # Connections
     ##################################################
     #connect square wave with multiply
     self.connect((self.analog_sig_source_x_0, 0), (self.blocks_multiply_xx_0, 0))
     # connect sine wave with multiply
     self.connect((self.analog_sig_source_x_0_0, 0), (self.blocks_multiply_xx_0, 1))
     #connect multiply output with audio sink
     self.connect((self.blocks_multiply_xx_0, 0), (self.audio_sink_0, 0))
     #connect multiply out to with null sink
     self.connect((self.blocks_multiply_xx_0, 0), (self.blocks_file_sink_0, 0))
     #connect multiply out to with null sink
     self.connect((self.blocks_multiply_xx_0, 0), (self.blocks_null_sink_0, 0))

     #if its non-zero send to file sink, otherwise send to null sink
       if  self.blocks_multiply_xx_0.out_sig != 0:
           self.disconnect(self.connect((self.blocks_multiply_xx_0, 0), (self.blocks_null_sink_0, 0)))
           self.connect((self.blocks_mgetinultiply_xx_0, 0), (self.blocks_file_sink_0, 0))
       else:
         self.disconnect(self.connect((self.blocks_multiply_xx_0, 0), (self.blocks_file_sink_0, 0)))
         self.connect((self.blocks_multiply_xx_0, 0), (self.blocks_null_sink_0 ,0))

def d():
my_top_block().run()
def n():
t0 = time.time()
time.sleep(5)
t1 = time.time()
t = t1 - t0
print 'time = ', t
def ram2disk():
#take the contents in written to Q: drive
#write these to a new file in C: drive
#execute this in a new thread occurring simultaneously as the write to ram

if name == 'main':
t = threading.Thread(name='non-daemon', target=n)
d = threading.Thread(name='daemon', target=d)
d.setDaemon(True)

d.start()
t.start()</pre>


texts: Timing syncronisation using polyphase decimation filter

texts:
This question is regarding application of polyphase decimation filter in symbol time correction. How does choosing a perticular "polyphase branch" in a polyphase decimation filter result in timing adjustment ?

Suppose I have a 256:1 decimation filter with a sampling frequency of $256\textrm{ kHz}$ with total of $256\cdot 8$ taps (or $256$ filter banks with $8$ tap each). These are linear phase FIR filters, each filter bank having a data rate of $1\textrm{ kHz}$. Consider following two scenario:

I start with polyphase filter bank index $1$ and carry on convolving input data and roll over polyphase bank index at $256$.
I start with polyphase filter bank index $2$ and roll over polyphase bank index at $256$.

How does this result in timing adjustment ?


texts: Window function for signal with 3-5 sine waves inside

texts: Should I apply some window before FFT if I want to find 5 frequencies in signal? Signal is pure 5 sinusoids + some noise generated during recording on analog cable.


texts: Is mean squared error or absolute error on two spectrograms of two different wav files enough to compare how similar they are?

texts:
I've been working on a genetic algorithm to figure out how to play a guitar chord just by listening to it. I have it all working now, but the system I am using to compare how 'similar' a guess is to the target seems less than robust. Here is the basic workflow of the algorithm:

I have a wav file and a corresponding spectrogram of the 'target' chord which I am trying to have the algorithm recreate
I create a population of 'chords'. I simulate how a person would place their fingers on a fretboard and then convert it to a more readable set of pitches.
Then, I run each chord through the simulation. I play what that chord would sound like through a MIDI library and, at the same time, record the sound to convert into a spectrogram. This is not ideal because then I have to wait for each chord to play in isolation, but to my knowledge, there is no way to just dream up a spectrogram, I have to just record it.
Then I assign each chord a 'fitness' value by doing either mean squared or absolute error pixel-wise between the target spectrogram and each chord's spectrogram. I also cut any silence before or after the wav file before turning them into a spectrogram and take the average across the time dimension for each spectrogram (so time is not a factor, and instead of being a plot of frequencies vs time, it is just a list of frequencies and their corresponding amplitudes)

However, I have noticed that this algorithm for finding the fitness of each organism doesn't work great. Even when I run the same chord shape through multiple times, I get a lot of fluctuation each time as to what fitness it receives. It should be relatively the same volume each time, so I can't imagine normalization will help, but I could try it. Any other ideas about how to get it more reliable? Are spectrograms not the way to go? Or maybe there is a smarter option than pixel-wise error? Thanks!


texts: When to consider double (64 bit) floating point for Audio

texts:
When synthesizing and processing audio on modern processors, when would one consider using something other than single precision (32 bit) floating point? Obviously, the audio coming in and going out to the real world is 16/24 bit, so I'm just talking about the precision of the signals (both the audio itself and things like filter coefficients) in the software.

Assume that:

the CPU/DSP has hardware floating point support for both single and double precision
The priority is quality audio, not high performance. For example, double precision would be considered if it offered better (perceptual) quality.


texts: What is the use of the Welch lower bound when k > 1?

texts:
The Welch lower bound on the maximum cross-correlation (inner product) of $m$ vectors, each of $n$ elements, is given by $c_{max} = f(m, n, k)$ (eg, see Wikipedia Welch bounds). It is always assumed that $m \gt n$ (otherwise the bound is 0).  Apparently the Welch bound has been widely applied in the restricted case that the vectors are real-valued and $k = 1.\;$ In that case, $c_{max} = \left[ (m-n)/(n(m-1)) \right]^{1/2},\;$ and is always real. The function has practical use in allowing us to compare the $c_{max}\;$ calculated from any given set of vectors with the theoretical minimum possible value.

In contrast, when $k > 1,\;$ then $c_{max}\;$ is frequently a complex quantity, even if the vectors are real-valued.  $\;c_{max}$ becomes real only if $m >> n\;,$ or more specifically, when $m \gt \binom{n+k-1}{k}.\;$
It is unclear to me whether $c_{max}$ has any use when $k \gt 1.$

$\quad$ Q1:  Does $k$ have any physical interpretation?

$\quad$ Q2:  Is the Welch formula with $k \gt 1\;$ useful in any area of physics, computational science,
$\quad\quad$ electrical engineering, or telecommunications? How do you interpret complex returned values?


texts: What is the relation between noise (AWGN) variance and number of antennas?

texts:
I'm working in two Matlab scripts trying to simulate a MIMO system using OFDM (downlink) and SC-FDE (uplink).

My question is: for different numbers of transmit and receive antennas, what is the relationship between the AWGN variance and the number of antennas?
Does it depend if it is uplink or downlink?

I'm using matricial equations like: Y = HX + N , with matrix dimensions: Y(NR x length) , H(NR x NT) , X(NT x length) , N(NR x length)


texts: Positive slope of BER with a greater frequency deviation?

texts:
I am trying to implement a frequency-modulated communication channel on MatLab. However, there is something I don't understand about the sensitivity of the frequency modulator (or the frequency deviation as my signal has an amplitude equal to 1). My message has a 500Hz frequency and is sampled at 20kHz. I added white gaussian noise to my bandpass signal (carrier of 6kHz).

Why does the BER of the demodulated message have a positive slope after 250Hz ?

NB : it is written [dB] for the BER on the graph but it is not, I forgot to remove it


texts: Can somebody explain SDF and MACE to me?

texts: I can't find any information on SDF and MACE save for a vague exert in Kinser's Image Operators books. Therefore I am dumbfounded. I can't understand what they do and what's their role in composite filters. Can somebody please explain them to me? Thanks.


texts:
I am working on the localization problem of an underwater vehicle. However the problem is still in very simple and it does not matter that it is about underwater.

How can I use the EKF when I have multiple measurements but with different sampling? Not all the measurements are available at the same time to perform the update step.

Can I perform the update for the measurements that I have at each time step? For example, I have a pressure sensor and a speed sensor with frequency $1$ Hz. But I have also another sensor (USBL) system with $0.01$ Hz.

What if I perform update step for pressure and speed sensor at each time step adjusting the measurement model to that case (like I have only these sensors) and once I receive the USBL measurement I augment the measurement vector including the USBL measurement?


texts: Beat detection and calculation of start time

texts:
Beat detection
I'm trying to create a beat detection function of a machine producing  in the language R. I hosted the audio file on dropbox.
It was recorded with 16 bit and 16kHz. I started with some filtering. The frequencies can be seen here. Most backgroundnoise is around 100Hz and I don't really have information above 4000Hz, so i used a bandpassfilter from 300Hz to 4000Hz. The filtered audiofile is also hosted on dropbox and the frequencies now look like this.
I then downsampled to 4000Hz to save data and computation time. Finally I used the hilbert function in combination with autocorrelation to get something like this.
With the peak I can calculate that the machine cycle is 70/min, which is correct.
This is the first time I'm doing something like this, so how can I improve the preprocessing? Should I downsample even more to save recources?
This method only looks at the loudness of the signal, disregards the frequency completly. Do you think thats good enough for something like this or can you recommend a method that also includes the frequency of the signal?
Calculation of start time
I also have to detect where the machine begins its production. I cut the signals into 1.5 second packets and labeled them. I used some simple metrics like mean and max of the loundness to cluster them like this. My idea was to check wheter its running or not, then use the packets of the running machine to calculate the cycle time. Then I can loop back the amount of the cycle time and check, if the machine is running at that point or not.
I don't really like that approach though. Do you have a better idea? Are there some better features than mean and max?


texts:
Im stuck on this one. I hope someone helps me!

Problem: I need to find Dopler shift with array of reflected chirps.

I have bursts of chirps which then i sampled, and calculated fft.
Im suposed to calculate Doppler shift using just spectrum. I see it's shifted but have no idea how to calculate the shift. I tried with boxcars in hope it will  be easier and i will notice something, but it's the same thing. Also i find some really useful article and also stuck understanding it.

  "Unless the target is moving at an extremely high speed relative to
  the speed of light, the Doppler shift will be small and very difficult
  to detect from one pulse. The solution to this problem is to transmit
  a burst waveform containing repeated pulses. Initially, for
  simplicity, boxcars were used as the pulse so that some insight into
  how the parameters of the pulse affect the Doppler shift might be
  gained."
  source


texts:
I have been trying to do the cross correlation of two signals, and I am expecting the output to have the same dimensions as the two signals have. However what I am seeing is that the cross correlation has twice the amount of points

The red here is the cross correlation of the green and blue. The two starting signals cover $1$ hour of time, sampled at $100Hz$, where as the cross correlation seems to represent $2$ hours, what I want to see is how correlated they are over the 1 hour. Can anyone explain this? (also apologies the graphs look like the joker)


texts:
I was reading this thread: How to get Instantaneous Magnitude for a Instantaneous Frequency From FFT?

I have basically the same question, but I'm curious about more detail about the Hann window kernel. ederwander said that once you know the frequency offset from the nearest bin then the instantaneous magnitude is a the Hann kernel of your offset times FFT length, times magnitude of the bin you're offset for.

Where the Hann window kernel $W(k) = 0.5 (\operatorname{sinc}(k(M/N)) / (1 - (k(M/N))^2)$

I'm just not sure why that kernel evaluated at offset*N gives the proper scaling factor for adjusting the magnitude. I couldn't even find any other information on the Hann kernel or how it is derived.

Thanks!!


texts: Cross Domain Equivalent to Nyquist Sampling Theorem?

texts:
In attempting to answer this question by @Oliver here: What characterizies 'causality' for a finite FFT? I have considered the minimum requirement to avoid time domain aliasing in the Discrete Fourier Transform, or more generally any application where the frequency domain is sampled. Similar to sampling in time at at least twice the highest frequency to represent the spectrum without the effects of aliasing, I suggest using a time duration that is at least twice the response time of the underlying continuous time signal to represent the continuous time domain signal (in the DFT) without the effects of time aliasing. Or when the time domain process is restricted to known causal processes, the time duration is at least as long as the response time.

This is the equivalent to Nyquist's Sampling Theorem in the frequency domain; ultimately "sampling in frequency" such that the duration of the time domain waveform is greater than twice its response time.

I understand that the same theory would apply, but given that Shannon in his paper provides the Nyquist theorem in the time domain specifically has made me curious if this property may go by other formally named theorems in other domains?

To illustrate this graphically consider the drawing by RBJ below except replace the frequency axis with the time axis.

Public Domain, https://commons.wikimedia.org/w/index.php?curid=1065579


texts: MATLAB : Proper estimation of weights and how to calculate MSE for QPSK signal for Constant Modulus Algorithm

texts:
I had asked conceptual Question about Constant Modulus Algorithm. I am implementing the simple steps of the algorithm, as I do not have the cma() built in module.

I am considering a Finite Impulse Response (FIR) system  whose true coefficients are $h = \begin{bmatrix}1 &0.45 &-0.2\end{bmatrix}$. The algorithm is presented briefly :

$u_n$ : Output of the FIR filter that is driven by Input signal, $s$ is a white Gaussian signal that drives an FIR process by the command

s=round(rand(1,N))*2-1;
u = filter(h,1,s);

We receive a noisy signal which is the output of the system corrupted by AWGN. Let the noise corrupted signal be

$$x_n = h^T s_n + \eta_n$$
Constructing an equalizer

$$y_n = w^T x_n$$

The cost function to be minimized by gradient descent is

$$J(w) = E\left[\left(\big\lvert y_n\big\rvert^2 -1\right)^2\right]$$

The weight update equation is given by

$$w_{n+1} = w_n - 2\mu e_n  y_n^T x_n$$

where error,

$$e_n = \left(\big\lvert y_n\big\rvert^2 -1\right)$$

Operator $T$ is the transpose and assuming real signals without any imaginary and complex parts.

Below is the code. The algorithm does not give proper result in terms of estimates of the weights.

Problem 1: The Graph for the Mean Square Error (in DB) for the estimated weights returned by the algorithm vs. Signal to Noise Ratio is giving an opposite trend i.e. instead of MSE decreasing with increasing SNR, I am getting the opposite when I re-run the program!! Below is the image of what I mean. The first figure is correct, but with the same code, I ran it again and got the second figure. Why is this & how can I prevent it?

I am not sure if this is due to the initialization of the variable $L$, filter order = number of delays = 2. Somethings are not straight for me which are for the FIR filter of the form : $u(t) = e(t) + 0.45e(t-1) - 0.2e(t-2)$, what is the filter order, smoothing length and the number of weights unknown to be estimated?
 Can somebody please help so that it works?

   clear all
clc

N = 256;
h = [1    0.45   -0.2];

R2 = 2;

mu =  1.0000e-009;
noisedB  =0;

L=2; % smoothing length L+1
ChL=1;  % length of the channel= ChL+1
EqD=round((L+ChL)/2);  %  channel equalization delay

i=sqrt(-1);

Ch=[1 0.45 -0.2]; %Channel
%Ch=[0.8+i*0.1 .9-i*0.2]; %complex channel
Ch=Ch/norm(Ch);% normalize

skip =1
for l=1:6
i          = 1;

 TxS=round(rand(1,N))*2-1;  % QPSK symbols are transmitted symbols
%TxS=TxS+sqrt(-1)*(round(rand(1,N))*2-1);
x=filter(Ch,1,TxS); %channel distortion

n=randn(1,N);   % additive white gaussian noise

n=n/norm(n)*10^(-noisedB/20)*norm(x);  % scale noise power
x1=x+n;  % received noisy signal

%estimation using CMA
K=N-L;   %% Discard initial samples for avoiding 0's and negative
X=zeros(L+1,K);  %each vector
for j=1:K
    X(:,j)=x1(j+L:-1:j).';  %y_n = w^T x_n
end

e=zeros(1,K);
w=zeros(L+1,1);
w(EqD)=1;    % initial condition
while i<=K

e(i)=abs(w.'*X(:,i))^2-R2;                  % initial error
   w=w-mu*2*e(i)*X(:,i)*X(:,i)'*w;  % update equalizer co-efficients
    cma_mse_h(l,i)    = sum((w'-h).^2)/3;
    est_w(i,:) = w;
   w(EqD)= 1;

            i = i+1;

end
noisedB = noisedB + 5;
end

for ii = 1:6

 Error(ii) = 10*log10(mean(cma_mse_h(ii,:)));
end

 plot([0:5:25], Error(1:6));
 grid on;
 xlabel ('SNR(dB)')
 ylabel('MSE_h')

Second Problem : The true FIR channel coefficients are not imaginary and have only real parts. But the estimated weights will have both real & imaginary if I work with real and complex representation. How can I properly calculate the weights & its MSE for this case?


texts: How to Extract Nonorthogonal PCA Principal Components

texts:
My problem is essentially a 'blind source separation' problem. I have 3 non-orthogonal sources (or basis functions) and $N$ random linear combinations (mixes) of said sources. My problem is to obtain the sources from the mixes.

Figure A shows the sources, B shows the mixes.

Approaches taken:

$\textbf{PCA}$ - I tried PCA on the mixes (fig.C), but the issue is that PCA will only give orthogonal bases, while my sources/bases are non-orthogonal. This issue with PCA is shown in figure D, where the data is clearly described by 2 non-orthogonal basis, but PCA (the solid lines) cant reconstruct them!
$\textbf{Factor Rotation}$ - I tried applying some solutions form Factor Analysis (not my forte). The promax rotation (matlab: nw = rotatefactors(cov,'method','promax') ) is shown in fig. D with the dashed lines. As far as I can tell, factor rotation works with the principal components, not the original matrix and thus I have no idea how it could reconstruct the right basis. I think this only works with factor matrices, not with generic ones...
$\textbf{ICA}$ - I tried overdetermined ICA with the fastICA algorithm (also not my forte, but I think I understand it better). I was hoping this would work since independent components (ICs) are non orthogonal. The solution is shown in fig.E. While the ICs are indeed nonorthogonal, they are NOT my original sources :-(

Any other potential tips or leads or solutions would be greatly appreciated.

More mathematical Detail

The 3 source signals are themselves linear combos of the 5 laguerre basis functions. Thus, the problem can be formulated as such:

I take l=5 laguerre basis functions of length m=30 to get B, an [l,m] matrix
I randomly mix then is my x=3 sources, S: S=B*M1, where M1 is a random [l,x]=[5,3] mixing matrix. S is displayed in fig.A
I mix S n=100 times to get my mixed signals, K (fig.B), where K=S*M2, a [x,n]=[3,100] randomly generated mixing matrix.
I have the original basis, L, & the mixed matrix, K. K=S*M2=B*M1*M2. My goal is to recover S.

NOTE: I understand that this may seem impossible since there are many possible S's to choose from (ie the factorizatio S*M2=B*M1*M2 is not unique), however, the true S, will dramatically decrease model complexity by having only x=3 rather than l=5 basis functions which can exactly reproduce all the kernels.

Thanks!


texts: Classifing audio signals to detect fault

texts:
The basic idea i'm working on is automatic fault detection using audio signals captured from Motor. I have like set of sample audio signal which are recorded when there is no fault and with fault at different rpm. While hearing the audio I can clearly identify and differentiate the audio signal. But how can I differentiate that technically or programtically. Basically im from computer science background , so im having a hard time to find the method of implementation whether it can be found by audio analysis itself or if I want to implement Machine learning on audio signals

Sample signals are here

1.3500rpm-clean
2.3500rpm-fault


texts: Indetermination in a complex fourier series

texts:
I determined the complex Fourier series of a sinusoidal signal and arrived at the following expression:
$$\sum_{n=\infty}^{\infty} \left[\frac{4e^{-j \frac{\pi}{2}n}}{\pi(1-n^2)}(e^{-j\pi n}+1)\right]e^{-j \frac{\pi}{2}tn} $$
However, we have a division by zero at $n = 1$ and $n = -1 $. Is there any way to determine $D_{-1}$ and $D_1$?


texts:
I have read through many papers related to the Mel-frequency cepstrum coefficient (MFCC): I have some confusion about the MFCC process:

Many of the papers say that I have to do a FT (Fourier Transform) on the log-spectrum to obtain cepstrum (as mentioned in Wikipedia : definition of power cepstrum)
Many say I have to do an IFT (Inverse Fourier Transform) on log-spectrum to obtain cepstrum.
Many refer to the DCT-II (Discrete cosine Transform-2) on the Mel-log-spectrum to obtain MFCC.

My question is, which one is right:

FT or IFT?
If the IFT is correct than why are they using DCT-II instead of IDCT-II, since the DCT-II is similar to DFT or FT?

(hence, the IDFT or IFT should be similar to IDCT-II)


texts:
There is already a succesful question How did they color old black and white movies? on this site. Very general and correspondingly in general answered. Let me please narrow it down to just the photograpy (i.e. no time dimension) and more specified starting point:

What is the algorithm to color the greyscale HD picture given the colors at various points on greyscale original?

Note: Yes, I am aware that the original colors can't be restored any more.


texts: Pose estimation of colored checkerboard cubes

texts:
I want to find the orientation of 2x2 checkerboard cubes, like these:

That is to say, I want to recover their positions and rotations within the image relative to the camera.

In my particular case, I care more about the rotation than the position (and the depth barely matters at all). I know the colors ahead of time, and can change them to be more convenient. I want to pick up the cubes, so a hand may be in the frame. The lighting may not be uniform, but the background should be simple (e.g. carpet, table, pavement).

I do have a partially working solution, which I've included below as an answer (it makes the question too specific to include it here). But I'm looking for suggestions on how to approach the problem and what techniques to apply.

Other things I've been considering:

Use a line detection algorithm to detect edges, and use those to seed hypotheses (instead of the center saddle point).
Do a cartesian-to-polar transform at the detected saddle points to turn the angled lines into straight lines (expensive?).
Score hypotheses by using them to unskew a face then dot product-ing that against a reference image of the face (but what about lightning changes?).
Random combinations of dilation and erosion throughout (invariably doesn't seem to help).


texts:
I try to extend my understanding of digital filter design, so I play around with some YouTube videos and the MicroModeler DSP with the following IIR design (based on a video):

I use the following parameters:

Zeros: 0.57 ± 0.78i
Poles: 0.51 ± 0.7i

How does the calculator generate the additional gain with value 0.90134? As I understand it so far the gain value is generated by the equation of the zeros, but I can´t calculate the value for it.


texts: Is R suitable for digital signal processing

texts:
While asking a question about representing large time series in R I was discouraged from using R for digital signal processing. I understand that R is geared towards statistics. However, a signal is just a set of measurements. R also deals with time series, but they are usually at different time scales. Nevertheless, e.g. when predicting price evolution it is a matter of understanding the model behind it and possibly predicting it. That's quite close to what I usually want in physics.

In "traditional" digital signal processing (DSP), it is more common to analyze the signal in terms of its frequency components and their interplay or coherence (which is something like frequency dependent covariance of the amplitudes of frequency components) between signals.

I'm interested in R mostly because of some packages on CRAN with advanced scientific routines (which aren't as developed in Python) like signal for general DSP, wavelets and biwavelet for wavelet transforms and timsac for bispectrum. These packages should make this more "traditional" approach possible. However, I get the feeling that time series are generally analyzed in a different way in R as is apparent from the time series task view on CRAN. Are these more statistical models perhaps unsuitable for analyzing non-stationary time series?

So I ask whether there is something that makes R worse for DSP than e.g. Python+SciPy (which I currently use the most) or MATLAB. Is it perhaps some underlying implementation that makes operations slower? Or is it just that it would be necessary to change one's approach to analyzing signals to a more statistical one when using R? And would such an approach be perhaps limited in some way in comparison to the "traditional" DSP approach?

Note:

I tried asking this question on SO, but I was directed to rather ask here.

I'm interested in objective answers with sources, not just opinions. E.g.

which internal representations makes it slower
what are the limitations of the statistical models used for time series analysis, e.g. how they cope with non-stationary signals, n-wave interaction, etc.


texts: Does the cv2.waitKey argument determine the frame-rate while capturing a video?

texts:
Kinda new to image processing, so pardon me if the question is too stupid.

Using the library cv2 in python to capture video from my webcam, which, so far as I can see as a programmer, is basically an array of images. But nothing specifies the frame rate of a video. Suppose I use the following code snippet to capture a video.

while True:
    capture, frame = cap.read()
    cv2.imshow('Being Captured. Hit Escape to Stop',frame)
    if cv2.waitKey(1)&0xFF==27:break

As I understand about the waitKey function, the camera captures a frame, waits for a millisecond, then captures the next one etc. So if a millisecond for each still image, does it mean the frame-rate is 1 KHz? Typical rates for most videos is 25-30 Hz, right? So can I adjust that just by changing the argument of waitKey in inverse proportion? Or, is my thought process wrong somewhere?


texts: Help with Audio Signal noise removing using FFT and IFFT in Matlab

texts:
I dont understand why my modified fft of a audio signal returns complex values after ifft.
I modified my signal by zeroing out all the unnecessary freq.

FFT: Fast fourier transfrom

IFFT: Inverse fast fourier transform

Heres the matlab code and graph

[a,fs]=audioread("Recording.m4a")

plot(a)

n=numel(a)/2  *//*since audio signal is stereo but i want only mono*

b=abs(fft(a))

y=fft(a)

c=y(:,1)     *//*to convert audio signal from stereo to mono*

plot(f,abs(c))

xlim([0 1000])

c(abs(c)<31)=0

c(round(n*1000/fs):end)=0

f=0:n-1

f=(f*fs/n)'

plot(f,abs(c))

xlim([0 1000])

z=ifft(c)

plot(z)


texts: Pitch detection to find music notes (from voice)

texts:
Hi signal processing community,

currently I am searching a topic for my master thesis in computer science and since I am a very good singer I want to do something with music.
I thought about a topic with pitch detection algorithms. But I am very new in this context and I want to inform me first, so I hope you can help me and answer me some of my questions.

I know that I have to work with FFT (Fast Fourier Transformation) and algorithms like AMDF or Yin. Can you give me a hint of a starting point (a paper or website).

I start with a signal which I got from my microphone, but what are the next steps.

I hope you can help me, thank you all.


texts: Cross-Correlation of two signals each with independent complex Gaussian noise

texts:
The nearest question I can find is here, but it does not have an answer. Say I want to correlate the following two noisy channels: $x_1 = s_1 + n_1$ and $x_2 = s_2 + n_2$, where $s_1$ and $s_2$ are arbitrary signals and $n_1$ and $n_2$ are independent complex Gaussian noise (identically distributed for starters). Further, let either $M \lt N$ or $M \ll N$, with $x_1 \in \mathbb{C}^M$ and $x_2 \in \mathbb{C}^N$. The cross-correlation of the two channels yields

$$
x_1 \star x_2 = s_1 \star s_2 + s_1 \star n_2 + n_1 \star s_2 + n_1 \star n_2.
$$

Perhaps I am searching with the wrong keywords, but I am having trouble finding information about $n_1 \star n_2$. What kind of distribution does it have? How much power does it contribute to the output?


texts: Why do butterworth filter sometimes make signals diverge?

texts: When using Butterworth filter in MATLAB, if I chose an order too high, sometimes the filtered signal would grow without bound to infinity, and if I decrease the order, the problem disappears, has anyone encountered the same problem, and why?


texts: Camera and laser scanner calibration final approach questions

texts:
A while ago I posted this question about camera and laser scanner calibration. I've been away from this project for a while and now I need to come back and get a final approach  to calibrate properly this system.

So having Cedron Dawg's as a good answer to get the laser plane and also having the method described in this article, I have the next approach (assuming I got already the camera intrinsic parameters and the distortion is corrected ):

Get the camera extrinsics (R|t) placing a chessboard in front of the
camera on the scanning area.
Get the next equations from this article:

Having the above functions obviously we need a third equation to match the number of equations and variables, so I though to use this answer approach to get the laser plane equation and with this I will be able to solve the three equation system for the world coordinates (X,Y,Z).

So assuming that for every camera frame I've got the laser pixel input image coordinates (x,y) I will be able to transform them to world coordinates (X,Y,Z) with the above equations.

is all of this correct? is there in this approach any mistake?

EDIT

I edit in order to clarify more what I'm trying to do. The next picture illustrates an example about what I'm trying to do:

The object will change in width (B) and height (A) uniformly ( assume laser,camera and target are stationary), so applying the laser I need to measure (for each laser point) height and width changes. So the aim is, for each camera frame, draw a calibrated laser profile of the object.

What would it be the best way to solve this?

Thanks in advance.


texts: Relationship between matrix rank and beamforming

texts: I always encounter the term matrix rank in papers related to beamforming. I am only familiar with the basics of beamforming (delay sum beamformer, basic capon). Can someone explain the significance of matrix rank in beamforming in a layman's term?


texts: Machine Learning in FAST corner detector

texts: This paper by Edward Rosten and Tom Drummond 2008 describes choosing interest point via machine learning approach using ID3 algorithm. I understand the concept of the ID3 they are using. But my question is how are they fitting this whole training and testing idea of machine learning, while are analysing real time video frame.


texts:
I'm working in a speech recognition project for the university. I'm using the TIMIT dataset and pylearn2 as machine learning library.
The project need to be simple, and it doesn't have to follow the state of art at all costs.

I've decided to use the RASTA-PLP technique to extract the features. More precisely I'm using the rastaplp.m function in Matlab. I want a sample for each phoneme, but when i execute the RASTA-PLP function on a single phoneme, i have multiple values for each phoneme. This is normal, considering that the RASTA-PLP function calculate his value over a certain number of temporal windows, and each phoneme has a certain lenght.
But i need for each phoneme, a single value. Actually, after having chosen the number of feature (14 in my case), i need an array (1X14), composed of one line and 14 cells.

How can a transform the result to that ? Is it ok just to calculate the average value over the different values of each feature ?


texts: Z-transform of an impulse signal in discrete time

texts:
am trying to compute the Z-transform of the following signal
\begin{equation*}
x\left[n\right]\:=\:\sum_{k=-\infty \:}^{\infty \:}\:\delta \:\left[n-k\right]
\end{equation*}
so I thought it would be
\begin{equation*}
X\left(z\right)\:=\:\sum _{n=-\infty \:}^{\infty }\left(\sum_{k=-\infty }^{\infty}\delta \:\left[n-k\right]\right)\:z^{-n}
\end{equation*}
exchanging the order of summation and including the z inside the parenthesis
\begin{equation*}
X\left(z\right)\:=\:\sum _{k=-\infty \:}^{\infty }\left(\sum_{n=-\infty }^{\infty}\delta \:\left[n-k\right]z^{-n}\right)\:
\end{equation*}
then we get
\begin{equation*}
X\left(z\right)\:=\:\sum _{k=-\infty \:}^{\infty }z^{-k}\:
\end{equation*}
which would be
\begin{equation*}
X\left(z\right)\:=\:\frac{1}{\left(1-z^{-1}\right)}
\end{equation*}

Is this correct?? can somebody help me please. Thanks.


texts:
So I need to program (using Java) some things that include generated sounds.
I am struggling to understand how to make it work because I require a rather short signal while still using a sampling rate supported by the average sound card.
For example a 100 uSec beep with a 200 uSec delay followed by another beep (initially supposed to be a click but a sine wave also works, would be grateful if you could explain why a sine wave makes it simpler).

My understanding is that I need a high sampling rate (number of samples per second), a frequency within hearing range (for example 1000) and a buffer to hold that data that will be fed to the audio stream provided by the Java libraries.
The issue is that with my limit of sampling rate (hardware limitations) I end up with an almost completely empty buffer because of the tiny intervals of sound it must hold to meet the time requirements.

Am I missing something? Is there another way?
Why does using a sine wave instead of a click make things easier?
Any help would be appreciated.


texts: Is there any transformation to exploit the sparsity of a Gaussian Wave?

texts: I am looking for a transformation in which the gaussian wave when transformed with a particular analysis function would make the energy contents be mostly present only in a short band of frequencies?


texts: Making the units of a analytically calculated PSD consistent with the units of an FFT

texts:
One thing that has always puzzled me is making the units of an analytically derived Power Spectrum Density (PSD) consistent with the units of an FFT.
Let's say we record the output of a frequency generator, which we just set to output a signal of the form
$$v(t) = V_{\rm{pk}} \cos (2 \pi f_0 t)$$
for some finite amount of time, $T$, and take the FFT of this  acquired signal. We might then expect to see something like this:

Now as I understand it the FFT is a Linear Amplitude Spectrum (LAS) which will have units of Volts, or, a Power Spectrum (PS) which will have units of Volts-squared, depending on what we choose do with the resultant FFT. But it is clear the FFT can only have units of volts because that is what the units of the time-transient signal are.
Now lets say I want to fit this FFT with some realistic function. This is straight forward as I have the transient signal and I can calculate, analytically a PSD from a time truncated Fourier transform:
$$\frac{1}{T}\left|\int_{0}^{T} v(t) e^{-2 \pi j f t} \ {\rm{d}}t \right|^{2}$$
which will be in units of Volts-squared per Hertz ($\rm{V^{2}/Hz}$).
So its clear that the units of the FFT and the units of my function I want to fit to it don't match up! Multiplying by some quantity in units of $\rm{Hz}$ fixes the problem for example either the bandwidth of the FFT, $\Delta f = 1/ T$ or the Effective Noise Bandwidth (ENBW).
How can i reconcile the units of my analytic expression and the units of my FFT?

The objective is to simply fit my FFT data, so I given that all this should only really affect the amplitude, I could of course allow all this to be absorbed by an arbitrary constant -- but it would be nicer to have a consistent result.


texts:
In a waveform designing i am getting this type of constellation(image_attached).
Can anyone please tell me why i am getting this type of constellation?
this constellation is after frequency correction.


texts: Wavelet transform and FFT using to extract feature power bands with EEG signals

texts:
I am using 5 channels [ fz , cz , c3 , c4 , pz] to detect drowsiness of driver
My First Question is, what is the right input to get feature power band ( Theta , alpha , gamma , beta ) to wavelet transform ? ( these 5 channels or 1 channel or what ? )
My Second Question is, Is it right to classify data based on theta only got from wavelet transform ?


texts: Meaning of Hop-Size in Filter Bank Interpretation of Short Term Fourier Transform

texts:
The Short Term Fourier Transform (STFT) is used to obtain the time-varying spectrum from a signal. There are two ways to understand it - the overlap add interpretation, and the filter bank interpretation.
The overlap-add method is where the signal is cut up into windowed segments in time and the DFT of each segment is taken. The window is advanced forward by R samples in time, where R is called the hop-size.

The filter bank interpretation is where the window is seen as a low-pass filter, and the segment is shifted in frequency by an amount equal to the time segment. The hop-size here is interpreted as a down-sampling in frequency.

I am unable to understand this interpretation and would like to know how R is interpreted as down-sampling of the signal in frequency. This would also clarify the constraints on choosing a hop-size R, to prevent aliasing in frequency (as having a hop-size R means down-sampling in frequency).


texts:
I have a data signal (column vector) with N values taken over at specific sample rate (dx) at a fixed frequency. After taking the Fourrier transform of this signal and representing it in the Fourier Domain (with fftshift), I need to eliminate the parts situated in the origin and in the positive side. The final signal I need is the one with the negative shift. Therefore, after filtering the unwanted parts I will need to perform an inverse Fourier transform to recover the signal in time domain and center it in the origin.
Finally, I do not understand how to filter the spectrum when there is no frequency variation and when I’ve tried specifying a sample frequency and then applying a matlab filter I didn’t get any results. Can anyone help me to understand how to make this work please?
PS: the spectrum is in the figure below


texts:
I am trying to solve the attached question. Please I would appreciate if some one can check my approach and tell me if I am correct and how can I proceed to find h0(t)

[![enter image description here][3]][3]


texts: Why does my simulation of the binary symmetric channel disagree with theory?

texts:
Could someone please answer me why the simulated performance (using matlab) of a Hamming code over a binary symmetric channel does not fit exacly with the theoretical?
Here is my code:

clear all

n=15;
k=11;
I=eye(k);

%define the parity matrix for the hamming code
P=[1 1 0 0;1 0 1 0;0 1 1 0;1 0 0 1;0 1 0 1;0 0 1 1;1 1 1 0;1 1 0 1;1 0 1 1;0 1 1 1;1 1 1 1];

%The generator matrix combination of parity and identity matrix
G=[I,P];

%Define parity check matrix
HI=eye(n-k);
H=[P',HI];

%verify that rows of H represent the dual space for the rows of G
%(verification = 0)
verification=rem(G*H',2);

UoB=1043162;

%convert UoB number to binary
ci=dec2bin(UoB);
ci=ci-('0');

%take the 11 rightmost digits
cw=ci(:,10:20);
cw=rem(cw*G,2);

%monte carlo simulation index
mc=10000;
index=1;
Pe_matrix=zeros(21,mc);
Pe_vector=zeros(1,21);

In the procedure demonstrated bellow random error patterns error=rand(1,15)<p; are generated for 21 different crossover probabilities.Consequently these patterns are subtracted from the original codeword c=xor(cw,error); in order to find the errored codeword (in any place where the error pattern has an 1 the codeword will flip its bit. After that, the probability of error for the current error patern is calculated Pe=(sum(xor(cw,c)))./k; by summing up the ones returned (the places where the errored codeword and the original differ) from the subtraction of the errored codeword and the original one,
and by dividing them by the total length of the codeword k=numel(c);.The Pe_matrix stores all the Pe values for each of the cross over probabilities and each of the mc iterations.

for i=1:mc

    for p=0:0.01:.2
        error=rand(1,15)<p;
        c=xor(cw,error);
        k=numel(c);
        Pe=(sum(xor(cw,c)))./k;
        Pe_matrix(index,i)=Pe;
        index=index+1;
    end
    index=1;

end

In order to obtain smoother results Monte Carlo simulations were performed. In each of the Monte Carlo iterations,each row of the Pe_matrix has the Pe values for the different crossover probabilities.Now by summing up the elements of each row and by dividing them with the total number of the Monte Carlo simulations we acquire a strong estimate of the BER
for the Binary symmetric channel.

for j=1:21
    Pe_vector(1,j)=sum(Pe_matrix(j,:))/mc;
end

%In the code demonstrated bellow the theoretical performance of the channel
%is calculated in order to compare it with the simulated.

for pe=0:0.01:.2
    index=index+1;
    Pblkerror=1-(1-pe)^15-15*pe*(1-pe)^14;
    Pd(1,index)=(3/15)*Pblkerror;
end

%finally by creating a legend one can compare the theoretical with the
%simulated performance by the depicted BER graphs.
semilogy(Pe_vector,'r-o')
hold on;
semilogy(Pd,'k-*')
hold on;
legend('simulation','theoretical performance');
ylabel('Probability of error,Pe');
xlabel('cross over probability');
title('BsC simulation');
grid on;


texts: How to remove noise from this image without causing the crack to 'shrink'?

texts:
I am currently embarking on an image processing project to detect cracks on surfaces using OpenCV and Visual Studio via C++. I have tried writing my program and at the moment, I am stuck because I have difficulties striking a balance between noise in an image and preventing the 'crack' from shrinking due to the dilate() operation. I am completely new to OpenCV and have only coded in C before, so I hope you will be patient and explain things slowly. :)

I am currently working on the photograph below, taken from Google Images, which has the parent URL: http://www.newindianexpress.com/states/kerala/article1330634.ece

I shall type my code below:

#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/photo/photo.hpp"
#include "highgui.h"
#include <cv.h>
#include <iostream>
#include <stdlib.h>
#include <stdio.h>

using namespace cv;
using namespace std;

int main( int argc, char** argv ){

    //Declarations
    double alpha; /**< Simple contrast control */ //unused
    int beta;  /**< Simple brightness control */ //unused
    float minVal;
    float maxVal;
    Mat bw_img;
    Mat contrast1;
    Mat gaussian_img;
    Mat filtered_img;
    Mat weighted_img;
    Mat adaptive_img;

    //Loads image from file and put image matrix into "original_photo"
    Mat original_img = imread ("crack.jpg", 1);

    //Convert color-to-gray image
    cvtColor( original_img, bw_img, CV_BGR2GRAY );
    //imshow("bw_img", bw_img);

    //Gaussian filtering
    int i = 7;
    GaussianBlur(bw_img,gaussian_img,Size(i,i),0,0);
    imshow("gaussian_img", gaussian_img);

    //HPF kernel
    Mat kernel = (Mat_<float>(3,3) <<
        0,  -1, 0,
        -1, 5, -1,
        0,  -1, 0);

    filter2D(gaussian_img, filtered_img,-1, kernel, Point(-1,-1),0);

    //Extract high-frequency features, e.g. edges  show up.
    addWeighted(filtered_img, 1.5, bw_img, -0.5, 0,weighted_img);
    imshow("weighted_img",weighted_img);

    //Perform adaptive thresholding
    adaptiveThreshold(weighted_img,adaptive_img,255,ADAPTIVE_THRESH_GAUSSIAN_C, THRESH_BINARY, 51, 20);
    imshow("adaptive_img",adaptive_img);

    //Construct structuring element
    Mat struc_elem = getStructuringElement(MORPH_RECT,Size(5,5));

    //Perform erosion first
    Mat eroded_img;
    erode(adaptive_img,eroded_img,struc_elem);
    imshow("eroded_img",eroded_img);

    //Perform dilation next
    Mat dilate_img;
    dilate(eroded_img,dilate_img,struc_elem);
    imshow("dilate_img",dilate_img);

    double otsu_thresh_val = threshold(dilate_img, dilate_img, 0, 255, CV_THRESH_BINARY | CV_THRESH_OTSU);
    double high_thresh_val  = otsu_thresh_val,lower_thresh_val = otsu_thresh_val * 0.5;

    Mat cannyOP;
    Canny( dilate_img, cannyOP, lower_thresh_val, high_thresh_val );
    imshow("cannyOP",cannyOP);

    //wait for long long time
    int c = waitKey(20000000000000);

    //if you get impatient, end program
    if (c == 27)
    {
    return 0;
    }
}

The images I got so far are rather noisy, with a lot of black specks after doing erosion followed by dilation. I shall post the photos later as I don't have enough reputation to post more links. Sorry.

I have tried to increase the blockSize and C values in adaptiveThreshold() function, from the current 51/20 to 101/50. The noises were mostly removed, but the crack feature got reduced in size. This image shall be posted later and referred to as Image 2.

I have also tried to maintain the adaptiveThreshold() figures for blockSize and C, and tried to perform a dilation followed by an erosion. The final image was less noisy, though there were still black specks, but the crack became slightly discontinuous. This image shall be posted later and referred to as Image 3.

May I ask, how should I remove noise via either adaptiveThreshold() function or the erode() and dilate() functions without causing the crack to become thinner or wider?

This is because if the crack is thinner, the crack will 'disappear' totally in the final output image.

Thank you! :)


texts: How to find a variance of sample sequence

texts:
I have a sequence such as

$$r[n] = y[n]v[n]$$

$y[n]$ and $v[n]$ are zero-mean and statistically independent.

I need to find a variance of $r[n]$ and show that it is white and equal to $\sigma ^2_y\sigma^2_v$.

I know that $y[n]$ has variance of $\sigma ^2_y$ and I also know that
$$Cor(v_n, v_{n+m}) = \sigma ^2_v\delta[m]$$

This is what I have so far :

Finding the variance of $r[n]$:
$$\sigma_r^2 = \mathbb{E}([r[n]-m_r)^2]=\mathbb{E}[r^2[n]] - 2\mathbb{E}[r[n]]m_r+m_r^2$$

Then

$$\sigma_r^2 = \mathbb{E}[r^2[n]] - 2m_r m_r+m_r^2 = \mathbb{E}[r^2[n]]-m_r^2 $$

I am stuck at this point. Could you suggest how do I proceed from here?

The process is WSS random process.


texts: Otsu binarization with value preservation

texts:
I'm trying to implement a modified version of the otsu binarization algorithm. I'm trying to binarize document images. But in the binarization procedure I want the object (in this case the text) to retain its original grayscale value while the background takes the value of 255; that is, white. I'm posting a sample image version that I found in a paper.

This is the original image:

This is the resultant image I want to obtain:

Could someone please tell me how to do it in Matlab?


texts: Channel Shortening for underwater acoustic channel communication

texts:
I'm new in this forum, hope that we can learn from each other well.  I need to ask about the channel shortening for underwater acoustic channel.

what are the common algorithms in that field? is there a common source code in MATLAB which is used to be compared with such new algorithms?

I don't see many papers in that field, why?

please help me,

thank you so much.


texts: How to Generate a vector map for histogram equalization

texts:
Assuming I have a gray scale image Im with $0:255$ levels of gray, I want to generate a vector map that will preform a histogram equalization on the image. I have read the wiki entry on the subject as well as the matlab function section and yet to understand. (Note that I need to mathematically generate the vector and extracting it from [J,T] = histeq(I); does not answer my problem)

I am assuming I am wrong but it is my understanding that I should use the following code:

mapVec=imhist(Im);
EqVecMap=mapVec;
for i = 1:size(mapVec)
    EqVecMap(i)=0;
    for j = 1:i
        EqVecMap(i)=EqVecMap(i)+mapVec(i);
    end
end
mapVec=EqVecMap;

What is the correct process for generating this vector?


texts: Treating the measured frequency response of a system as a DFT

texts:
I have frequency response data of an analog system that I would like to turn into a digital filter in Matlab.  Is it kosher to think of my frequency samples as constituting a DFT of the system impulse response, even though there was never a time-domain signal in the first place?  If I do an ifft on the frequency response data, can this be treated as the impulse response of the system?  If so, what is the sample rate?  Do I assume it is twice the highest measured frequency?

Is there a way to build a digital filter from impulse response data without knowing the sample rate?


texts:
I have data of different SPL for different 1/3 center frequencies.
Is it correct to sum the pressures at a mean square basis of those SPL values (and then convert them back to SPL) to obtain the OASPL (Overall Sound Pressure Level)?
Furthermore, how to compute a SEL (Sound Exposure Level) from different OASPL over a time period? Just acoustically sum those OASPL values and divide the sum by the number of time steps?
Many thanks in advance!


texts: Probability of Detection and pulse-pulse decorrelation time

texts:
I'm trying to analyze the probability of detection ($P_d$) for a Swerling II target. I know that the Swerling II target assumes that the Radar Cross Section (RCS) is an exponential probability density function, and that the RCS fluctuates from pulse-to-pulse.

My question is: Do the Swerling models assume perfect decorrelation from pulse-pulse? As in, do they assume that the RCS on each pulse is purely independent of other pulses, or do they account for some correlation?

The target I'm trying to analyze does not perfectly decorrelate between pulses. So, the RCS on pulse #2 will be somewhat correlated to that on pulse #1.


texts: Is it customary to correct for the gain of a window?

texts:
Consider how the Hanning window is defined:

0.5 - 0.5 * cos(n*2*Pi/(N-1))

By this definition, it has a gain of 0.5, which is simply the average value of the coefficients. By contrast, Flattop windows, as defined, have unity gain, presumably by design.

It would seem appropriate to scale the Hanning window by a factor of 2, but I have never seen this discussed anywhere. It would seem that all windows should be scaled for unity gain.

In practice, are windows typically corrected for their gain?  If not, why not?

EDIT:

Since nobody has given an answer, I'll elaborate a bit.

It is quite easy to find papers that report the gain of the more common windows. But nowhere have I seen anyone refer to correcting the gain before using it for spectral analysis. Maybe I have always missed that statement, or everyone assumes gain correction to be an obvious requirement.

It seems like common sense to set the gain of a window to unity so that signal's energy level is preserved. Furthermore, how can one compare the various windows for amplitude accuracy if one has 0 dB gain, as a flattop does, and the other has nearly 10 dB loss, as the Gauss does.

Windows are also widely used for FIR filter design. In this application, it should be clear that the signal to be windowed, a sinc pulse, has most of its energy in the center of the window. Consequently, the window does little to reduce the sinc pulse's total energy. Thus, when used for filter design, we don't want unity gain, but rather unity peak amplitude, as most windows have, except the flattops. Something other than unity peak amplitude would affect the gain of the resulting FIR filter.


texts:
Assume a monostatic radar is transmitting a pulsed waveform. The pulse has been modulated by a Barker code with$\ N$ chips or bits. The waveform chip rate is $\tau_c $, the radar receiver's detection bandwidth is $\beta_c = 1/\tau_c$, and the waveform's pulsewidth is $\ N*\tau_c$.

I've heard conflicting opinions about a signal processing gain,$\ P_g$, due to compressing the pulse that is proportional to the waveform's time-bandwidth (or pulse compression ratio (PCR)) product (or equivalently the product of the uncompressed pulsewidth $\tau_p$ and the compressed bandwidth $\beta_c$) of the pulse: $\ P_g = \tau_p*\beta_c=N\tau_c*\beta_c$

Assuming that you start with the power form of the radar range equation:

Is it appropriate to add a$\ P_g$ term to the power form of the SNR radar range equation?

(Pulse compression allows a radar to transmit a longer pulse to get more energy on the target while simultaneously achieving the range resolution of a much shorter pulse.)

FWIW, I have the same question about the energy form of the radar range SNR equation:

My question is similar to this older question but more specific


texts: In the book digital Signal Processing by Proakis 4th Edition, at page 593 a table is drawn about some second order modules for discrete time systems. In that table there Transposed Form II drawn But when I searched the net, it was drawn completely opposite to that drawn in the book. On the internet, the input and output were placed on right and left respectively. I want to ask that which form is correct?


texts: How to calculate transfer rate and needed bandwidth

texts:
I have two questions:

1- If I want to transfer a photo via bluetooth, how should I find how long it takes to transfer? I mean the formula between the image size, bluetooth bandwidth and time.

2- If I want to transfer video in real time, how should I calculate the bandwidth is needed to transfer the video via bluetooth/WiFi? What's the formula between video resolution(VGA,2K,4K,etc) and transferring bandwidth if we want to do it in real-time?

For example: If I have 1920*1080*30(fps) video, I think I should multiply it by 3 or 4(respect to colored image type) then I have about 180~240 million bytes to transfer per second. So it means I need about 200MBps bit-rate! And I only should use 802.11n protocol for that!

Did I calculate the things right?

EDIT: I found this calculator that I am not sure but I think it's the thing I was looking for (https://www.cctvcalculator.net/en/calculations/bandwidth-calculator/). what do you think about this? Why does it calculate local bandwidth bigger than network bandwidth? it's strange for me!


texts: Laplace transform : integral vs poles & zeros

texts:
If Laplace transform is expressed as :
$$\int_{-\infty}^{+\infty} h(t)e^{-st}dt $$
with :
$$s = \sigma + j\omega$$
and $h(t)$ an impulse response expressed as :
$$h(t) = Ae^{-\sigma_0t}\cos(\omega_0t+\phi) = e^{-\sigma_0t}\cos(\omega_0t)$$ ($A=1$ and $\phi = 0$ for simplification, $h(t)=0$ if $t<0$)
Then, each vertical line (parallel to the imaginary axis) in the $s$ plane corresponds to the Fourier transform of $f(t) = h(t)e^{-\sigma t}$ for a fixed $\sigma$.
For $\sigma = -\sigma_0$, the decaying exponential of $h(t)$ is canceled and we get the Fourier transform* of $h(t) = \cos(\omega_0t)$, that is : diracs at $\omega_0$ and $-\omega_0$ (not accurate, see (*) just below), hence two poles : $-\sigma_0 + j\omega_0$ and $-\sigma_0 - j\omega_0$ as in the following picture (illustration only, poles not located correctly) :

Indeed, we can understand that :
(*)Please, note that the following is not accurate : since $h(t) = 0$ if $t<0$, we should use the unilateral Laplace transform, not bilateral ! So here we would get the unilateral Fourier transform of a sinusoid, not the bilateral (with diracs only) one ! To see what this would be, please see the link given at the end of the accepted answer
$$\int_{-\infty}^{+\infty} h(t)e^{-j\omega t}dt $$
$$= \int_{-\infty}^{+\infty} \cos(\omega_0t)e^{-j\omega t}dt$$
$$= \int_{-\infty}^{+\infty} \frac{e^{j\omega_0t}-e^{-j\omega_0t}}{2}e^{-j\omega t}dt$$
$$= \frac{1}{2}\int_{-\infty}^{+\infty} e^{j(\omega_0-\omega)t}-e^{-j(\omega_0+\omega)t}dt$$
If $\omega = \omega_0$ or $-\omega_0$, then the integral would blow up due to the $$\int_{-\infty}^{+\infty} e^0dt $$ member, hence the poles in the s plane.
So as shown in ch.32, p.24 of The Scientist and Engineer's Guide to DSP (see figures below), with Laplace transform we multiply $h(t)$ with $e^{-st}$ = $e^{-\sigma}e^{-j\omega}$, that is we multiply $h(t)$ with sinusoids that are either :

(a) Exponentially decaying ($\sigma$ > 0)
(b) Stable ($\sigma = 0$)
(c) Exponentially growing slower than our impulse response decay ($ -\sigma_0 < \sigma < 0$)
(d) Exponentially growing, compensating our impulse response decay ($\sigma = -\sigma_0$) : OK, as studied above.
(e) Exponentially growing quicker ($\sigma < - \sigma_0$ and $\sigma < 0$)

(letters correspond to pairs of points in the s plane shown in figures below, always at a fixed $\omega$ or $-\omega$ value)

I understand case d : since we cancel the exponential part, we get only the (unilateral !!) Fourier transform of a sinusoid. That is : infinite at $\omega_0$ and $-\omega_0$ hence the poles (though I don't know why we have a continuous function of omega with infinite values at $\omega_0$ and $-\omega_0$ instead of diracs as in the original Fourier transform of a sinusoid -> Because we use unilateral Laplace hence Fourier, see end of accepted answer !).
Case a, c and e are intuitive. In case a, we multiply $h(t)$ with a decaying exponential. The integral will be some finite complex value (for all values of $\sigma > 0$. In case c, we multiply by an exponential growing slower than the decaying exponential of $h(t)$, hence some finite complex value for the integral (for all values of $-\sigma_0 < \sigma < 0$). In case e, we multiply the $h(t)$ by an exponential that grows quicker than exponential of $h(t)$ decays : hence integral does not converge (for all values of $\sigma < -\sigma_0$).
But for case b, I can't get the intuition of why the integral would be zero as shown with the area under the curve (red in the above figures) ? In other words, I understand the vertical line in the s plane at $\sigma = -\sigma_0$, it is the Fourier transform of $h(t)e^{-\sigma_0 t}$ so it is Fourier transform of $h(t)$ once its exponential component is removed, hence 2 poles due to sinusoid. We get poles whenever $e^{-st}$ is identical (compensates) to the impulse response. But what would cause Fourier transform of $h(t)e^{-\sigma t}$ to be 0 at some $\omega$ ? For which $h(t)$ and how it would impact the area under the curve (integral) ?


texts: Fourier Series of Aperiodic convolution of periodic functions

texts:
we were given the following classic exercise:

  Given two periodic signals $x(t), y(t)$ with fundamental period $T$
  with Fourier series coefficients $c_m^x, c_m^y$ respectively, find the
  Fourier coefficients of the signal $z(t) = x(t) * y(t)$ with relation
  to $T, c_m^x, c_m^y$.

Now, this can easily be solved when the aforementioned convolution is the circular convolution (integral over a period only).

However, in class our professor noted that it can be solved even when we have an aperiodic convolution (that is, convolution as an integral from $-\infty$ to $+\infty$). We argued that, in this case, that infinite integral doesn't converge, and he responded that, even though the convolution integral doesn't converge (i.e. might be infinite), the Fourier Series coefficients are still finite and can be calculated!!

Is this true? If yes, then is the relation the usual one: $c_m^z=Tc_m^xc_m^y$ or another and how do you prove that? If not, why? Austere mathematical proofs would be appreciated.


texts: How to remove certain background from image?

texts:
This was the image obtained after Gabor filtering... Is there some way of removing the lines inside the image except the bright white rounded defect.

Did try another orientation of gabor filter go the following result:

After thresholding :Which is ok enough but not great

After Median filtering:

Thanks in advance


texts: How can I change parameters into a text string in MATLAB?

texts:
I am trying to compute the Morlet wavelet transform (cmor) of a signal, but I have to do it several times with different parameters Fb and Fc. Is there a way to automate it? For example, I have to do the transform with: cmor1-1, cmor1-2, ..., cmor1-10.
I created a vector Fc=1:10, but I'm not able to replace Fc parameter into cmor1-Fc.
I tried this way:

morlet='cmor1-Fc'
for i=1:10
morlet(1,8:9)=i;
end

but it doesn't work (there appear squares instead of the number, and it would properly work only for 2-digit numbers.
Thank you


texts: Random Signals - statistical properties are time dependant?

texts:
I'm taking a course on DSP and we're being introduced to the random signals, in particular continuous time and discrete time random signals.

We're told that if we repeat a single random experiment at each point in time, the random variable associated with each random experiment when collected together will form the random signal.
We are then told that the random signal/process has statistical properties such as mean that are time-dependant.

Now that makes no sense. Surely, if you're doing the same experiment, the mean of the entire random signal will be the same as the mean of the random variable for one experiment? How will time affect the random signal?