texts: Estimating Signal-To-Noise ratio from discrete time series

texts:
A hidden variable $x(t)$ (in my case, neuronal signal averaged over a brain area) exhibits dynamics on multiple timescales from milliseconds to seconds. An observable $y(t)$ is produced by convolving $x(t)$ with some kernel $k(t)$ which is in general unknown (in my case, $y(t)$ is a fluorescent signal of indicator molecules located inside some neurons, and $k(t)$ is the fluorescence response function to changes of neuronal calcium concentration, which is $x(t)$). Finally, $y(t)$ may acquire some instrumental noise $\nu(t)$, which for the purposes of this question will be assumed to be time-independent and normally distributed with unknown variance. $y(t)$ is sampled with 50ms intervals. About $k(t)$ it is known that it is causal ($k(t) = 0, \forall t<0$), that it ramps up quickly and decays slowly, and that its timescale (e.g. std) is somewhere between 200-600ms. $k(t)$ is typically approximated by an exponential or log-gaussian distribution.
My goal is to estimate SNR due to $\nu(t)$. I don't really care in what units is SNR measured, I just want to be able to discriminate noisy experiments from less noisy ones. Since white noise has peaked autocorrelation, I am considering using time-shifted autocorrelation as measure of SNR, namely
$$AC_1 = CORR(y(t), y(t+\Delta t))$$
where $\Delta t$ is one timestep. I have tested this metric on simulated data, where $x(t)$ and $\nu(t)$ were both gaussian random variables, and $k(t)$ was an exponential function, and found that $AC_1$ on average grows monotonously with increasing SNR, which gives me hope that it is a good metric for my task.
Questions:

Are there better metrics than $AC_1$? I have considered computing ratio of PSD for fast vs slow parts of the spectrum, but I am not quite sure what threshold to use, or whether it is actually better for my task
This whole analysis rests on the assumption that the autocorrelation of the original signal is "sharper" than that of the convolved one. In my simulations, if $x(t)$ was a very fastly changing function, $AC_1(y)$ would still be around 0.95, and for slower $x(t)$ it would be even larger. Is it possible (and realistic) that the convolved signal would actually be faster than the original? In that case, my analysis would wrongly mark fast signals as noisy, which is undesirable.


texts: Is it feasible to try and classify an enclosed space as being a. Very crowded b. Less crowded & c. Sparsely crowded using image processing?

texts:
The idea is to get three classifications for crowd: High, low and medium based on an image captured inside a train compartment. Something like this:

This gets classified in the 'High' category

http://www.spiceflair.com/wp-content/uploads/2012/05/14a.jpg
This perhaps in the 'Medium' & so on.


texts: From IQ signal to FM modulated carrier, how its done?

texts:
I am trying to comprehend GMSK modulation/demodulation. IQ samples are being phase keying is such a way that after multiplying them by complex carrier and sum both I and Q together, there is FM modulation. I don't understand how its being done. I guess it must be trigonometry but can someone explain me how phase changes in IQ result in frequency changes in HF signal?
Thanks in advance


texts: Is there a better Fourier Transform-based algorithm to use in Image recognition

texts:
I have been trying to use the Fourier transform in recognizing images of the same size 200x400.

I have tried many different ways to do that such as:

Performing the Fourier on the full image and then try to match it with the ones in the database. the accuracy rate was very bad.
Performing the Fourier on the full image and then take only the
sub-matrix contained in centre of the image and then match it with
the ones in the database. the accuracy rate was also very bad.
Performing the Fourier on the full image and then take only the
horizontal and vertical vector that intersect in the centre of the
image, and then match it with the ones in the database. also, the
accuracy rate was very bad.

Is there any better algorithm that might help me solve this problem?

Thank you very much


texts:
Suppose you have two designs of IIR filters with different specifications which are meant to be implemented on a DSP. On what parameters will one compare the two filters so that we can choose which filter is better (like Speed, Memory or Efficiency)? And how do these parameters vary filter to filter?

For example-
I have 2 filters, A and B:
Filter A is a 4th order Low Pass Elliptic filter and filter B is a 70th order Butterworth Low Pass Filter and both can do the required job satisfactorily.
Which filter should be implmented on the DSP and why? (Not taking into account the difficulties one might face in implementing a different algorithm)


texts: Convolution theorem with a kernel smaller than the image

texts: Suppose that I need to apply a convolution filter to an image. I want to do it by using the convolution theorem so I compute the kernel for the size of the input image to later calculate fft's and multiply. But I want the fft calculation to be less time consuming. The kernel loses its magnitude the further from its center, as usual. If I want to crop it, say, have a smaller kernel with only significant values, can I for instance compute its fft and later somehow extend the resulting frequency domain image of the kernel to fit the size of the input image, to later multiply them point-wise?


texts: Is concept of "bit" in computer programming similar to the concept of "bit" in information theory?

texts:
until today I knew that one bit is a variable, or a space in memory that can hold a value of either One (high) or Zero (low). This is the concept I learned from studying computer programming, microprocessor or DATA bus etc.

But after starting the course on information theory, I found out that bit is expressed as the information content of a symbol in message. This is calculated taking the logarithm (base 2) of the inverse of the probability of occurrence of the symbol.

Are these two concepts same ? On one hand one bit is a variable that can store either zero or one. On the other hand, one bit is the uncertainty associated with one of two symbols with probability of occurrence of 0.5. So, does 1 bit in computer programming or ASCII code mean 1 bit in information content of source or information theory?

A little edit:
here is one thing I am finding trouble understanding this topic. See, in data transfer of English alphabets, if we use ASCII code, we basically represent each symbol with 8 bits. Suppose that's 00000000 for a, 00000001 for b etc. So we are essentially allocating 8 quantization levels for each symbol.

But when the information theory comes into play, we take the probability of each symbol into account. 'E' has the highest frequency, where 'Z' has the lowest. So average information content comes down to 3 or 4 bits, right ?

My book says, 'Entropy or average information content is the minimum average number of bits required to represent each sample without distortion'.So, in this case, for efficient data transfer, are we creating maximum four quantization levels for each symbols? Because, on an average they carry information worth 4 bits. If that's so, isn't bit in information theory the same as the one in computer programming, data transfer or ASCII code etc ?

You probably get that I am clearly a noob here :p


texts:
I am trying to find the fundamental frequency of a low-frequency signal. I need an estimate that is precise to .01 Hz based on only a few cycles, so I'm trying to code up a fft in Python. The signal itself will be almost a pure tone, so I'm testing it out on a simple simulated wave.

I'm not seeing the results I expect, which may just be a dumb coding mistake. I'm new enough to DSP, though, that it's also possible I'm fundamentally misunderstanding what the transform is doing. Any help would be very much appreciated!

from math import *
import numpy as np
import matplotlib.pyplot as plt

I start by setting some system parameters

fs = 10000 # sampling rate (Hz)
T = 0.1 # length of collection (s)
assert fs*T % 1 == 0 # make sure we have an integer number of samples
windowlength = int(fs*T) # total number of samples

Then I create a sampled cosine wave for testing. Its period isn't an integer multiple of my window, so I'll need to zero-pad

f0 = 20.1 # fundamental frequency in Hz
checkwave = [0]*windowlength
for x in range(windowlength):
    checkwave[x] = cos(2*pi*f0*x/fs)

Plotting my wave in the time domain looks fine

t = np.linspace(0,T,int(fs*T)) # time domain x axis
plt.figure()
plt.plot(t,checkwave)

I then carry out the zero-padding and take the dft with numpy's rfft module

zp = 10 # factor to use in zero padding
dft = np.fft.rfft(checkwave,zp*len(checkwave))
fmaxindex = np.abs(dft).argmax() # find the peak frequency's index
fmax = fs*fmaxindex/(zp*len(checkwave)) # calculate the peak frequency
print(fmaxindex) # Higher than I expected
print(fmax)

But my plot in the frequency domain is off. I thought I should be looking at a sampled sinc wave centered on 20.1Hz, and since I'm sampling well above the Nyquist limit the centering should be perfect, right? But, if that were the case, 20Hz should be the largest bin (translates to bin 20). Instead, the largest bin is 21Hz. Plotting, we can confirm that's the case.

f = np.linspace(0,fs/2,int(fs*T*zp/2)+1)# frequency domain x axis
freqspan = 4 # number of samples to plot around the peak, so we get some detail
plt.figure()
plt.plot(f[fmaxindex-freqspan:fmaxindex+freqspan+1],abs(dft[fmaxindex-freqspan:fmaxindex+freqspan+1]))

My first thought was I made an off-by-one indexing error, but when I reproduce my last two blocks of code with an extra factor of ten on the zero padding, I see that things get even worse.

zp = 100 # factor to use in zero padding
dft = np.fft.rfft(checkwave,zp*len(checkwave))
fmaxindex = np.abs(dft).argmax() # find the peak frequency's index
fmax = fs*fmaxindex/(zp*len(checkwave)) # calculate the peak frequency
print(fmaxindex) # Higher than I expected
print(fmax)

Plotting in the frequency domain:

f = np.linspace(0,fs/2,int(fs*T*zp/2)+1)# frequency domain x axis
freqspan = 6 # number of samples to plot around the peak, so we get some detail
plt.figure()
plt.plot(f[fmaxindex-freqspan:fmaxindex+freqspan+1],abs(dft[fmaxindex-freqspan:fmaxindex+freqspan+1]))

Now the peak frequency is 20.8 (bin 208), about 6 bins away from where the sinc should be putting it.

Thanks in advance for any thoughts.


texts: Is a real Grayscale image always 8 bit per pixel indexed?

texts:
I tested several Grayscale images downloaded from the Internet. Some of them are 8bpp indexed while others are 24bpp (true colored). Looks like different Grayscale images have different depths.

For instance, this online application converts a color image to Grayscale. But, when I check its bit depth, I see that bit depth is not changed.

I think, Grayscales with 24bpp are pseudo-grayscales. That is, they are the average of R, G and B components. I think, real Grayscales are always 8 bit per pixel indexed.

Is my notion correct?

If not, please, provide me some links that I can study.


texts: Technical difference between bent pipe and regenerative transponders

texts:
I am preparing a proposal for satellite transponders and need to choose transponder architecture.  Based on primary research I realized that the regenerative transponders have a better characteristic.
I am looking for a reference to extract the transponder characteristic(for example BER and SNR relationship).
Does anybody have any idea?
thanks a lot


texts: Circular Time Delay in signal, effect on phase spectra of DFT

texts:
If you circularly shift a signal x[n]= [-3 -2 -1 0 1 2 3 2 1 0 -1 -2 ] to the right by M=1 compared to M=2, the phase spectra has a lot many zeros when M=2 compared to M=1.

here, I am talking about Circular Time shifting.

I know that if x[n] is even, then DFT is purely real and hence phase is zero. However, why are there way many more zeros when M=2 compared to M=1?

Thanks


texts: Remember that Christopher Walken skit on Saturday Night live where he told the band “it needs more cowbell” https://en.wikipedia.org/wiki/More_cowbell ?  This site automated it http://www.morecowbell.dj/ using echonest http://developer.echonest.com/ .  That’s all fine.  Question is: “What is the actual algorithm behind the beat detection and cowbell insertion?”


texts: Is the output of function of two ergodic processes ergodic?

texts: Let $\{\xi_k\}_{k\in \mathbf{Z}}$ and $\{\epsilon_k\}_{k\in \mathbf{Z}}$ be two independent zero-mean Gaussian processes (i.i.d.). Is the output of the function $f$ such that $y = f(\dots,\xi_{k-1},\xi_{k},\xi_{k+1},\dots,\epsilon_{k-1},\epsilon_k,\epsilon_{k+1},\dots)$ is an ergodic process? I found similar claims in the book ("Stochastic Processes and Long Range Dependence"). Can anybody show me some references or guide me through mathematical language that is easy to understand? Thank you!


texts:
I am working on IIR filter design, I got filter coefficients, now I want to truncate the coefficients so as to perform scaling.
let 0.2145678340000 be one of my coefficients. I want to truncate to 0.21456
I used below code to do that, but it wont work.
x_tr = fix(x*2^5)/2^5

can someone help me
thanks


texts: Standard deviation of Morlet wavelet's Fourier transform

texts:
MATLAB has the complex Morlet Wavelet in the following form:

$$\psi(t) = \frac{1}{\sqrt{\pi f_b}}e^{\frac{-t^2}{f_b}}e^{j2\pi f_ct}$$

I arrived at its Fourier transform as shown below (another arrived at the same: Fourier Transform of Morlet wavelet Function?).

$$X(f)=\mathcal F(x)(f)=e^{-\pi^2 f_b(f-f_c)^2}$$

I sought a relationship between $f_b$ and the standard deviation of the Gaussian form of $X(f)$. Is it right? Is there any citable published literature (book or journal article) arriving at this relationship?

Rearranging this into the familiar Gaussian form,
$$X(f) = e^{\frac{-(f-f_c)^2}{2\left({\frac{1}{\pi\sqrt{2f_b} })}\right)^2}}$$

The standard deviation $(\sigma_f)$ of this frequency domain Gaussian is then,
$$\sigma_f = \frac{1}{\pi\sqrt{2f_b}}$$
Rearranging it,
$$f_b = \frac{1}{2{(\pi\sigma_f)}^2}$$


texts:
This question is for people working in multi rate signal processing. Using compressive sensing techniques, we can still recover a signal by sampling below nyquist rate(by decimating the signal). How is that possible, can you explain?

There are some answers available on the web , but I want to hear it from SE group.

phani tej


texts: Direct Form I Filter from Direct Form II


texts: Applying a Liftering function to an audio signal

texts: Is it necessary to perform the liftering function before applying MFCC to an audio signal?


texts: What is the best approach to classify an image containing either text or image?

texts:
I am doing a research about image-based User Interface (UI) Analysis, and while UI image-based analysis is not more complicated than normal image analysis, I was able to do a segmentation process to separate text, icon and picture in the UI.

However, my next step is to perform content classification on the segmented groups.  The segmentation contains:

text
image or icon

So: what could be the best approach to classify an image containing either text or image?


texts: If two signals have different distributions that you know and can model, can you use that to separate them?

texts:
For example, if I have a mixed signal composed of signal A added to signal B, and I know that the histogram of signal A is non-Gaussian in the time domain, and that the histogram of signal B is Gaussian in the time domain, is there a way to isolate either A or B using this information?
I understand that I can use MLE or some similar method to fit a Gaussian to signal B, but I'm not sure where I'd go from there. My initial idea was to create a Gaussian filter that matches the distribution of signal B, and use that to separate B from A+B, but I'm not sure how I would implement that in practice, or if that even makes sense.
Any insights would be appreciated! Thanks


texts: How to add the noise covariance matrix of my measurements to tmy 1D kalman filter?

texts:
So I am new to kalman filters. This is my first experience with it.

I have some measurements but sometimes due to external factors those measurements become totally wrong (hence the peak when you run my matlab code). If I understand the principles of kalman filtering correctly, I can compensate this incorrect measurements by setting the measurement covariance matrix R as a variable, where R represents the accuracy of the measurement devide. But so far I am not sure about where I should add or how I should include this measurement covariance matrix in my measurements.

Could anybody help me out?

I have tried as follows: https://pastebin.com/KQEx5T10 But am totally not sure about my results.To me it looks like both of the outputs of my kalman filter are too far away from the actual measurements once the peak happened... I would have expected the output curves to be between the "toothshaped" curve and not above it.


texts: To characterize a signal it's a common use to determine bandwidth at -3dB from the peak (half power). Also, filter's cutoff frequency is pointed at -3dB. Why we use this rule?


texts: Find smoothed first derivative from signal with noisy slope

texts:
How can I get filtered first derivative from a noisy signal that has slowly changing slope in form of y=kx+b? k can slowly change in time, and I want to estimate its value.

I have tried 3 different approaches:

Take derivative as dx(i) = (x(i)-x(i-99))/100
Smooth with sliding mean window = 100, then take derivative as dx(i) = (x(i)-x(i-99))/100
Simple IIF filters (e.g. y(i) = 0.99*y(i-1) + 0.01*x(i), then take derivative as dx(i) = y(i)-y(i-1) and again filter with similar IIR, e.g. dy(i) = 0.95*dx(i-1) + 0.05*dx(i)

Problems:

Least-squares, regression and FIR filters (except rectangular window) have high computational cost, since I have to translate it to micro-controller with no DSP. That is why I can use only rectangular windows and IIR filters (they have low order).
If I find first derivative first, then smooth, it will be very noisy. So, I should smooth the original signal first, then find derivative from smoothed signal (and perhaps smooth the derivative again!).
I should play with filter parameters manually and it is hard to understand the frequency response of the whole system.

Question:

Maybe there is a single special (optimal?) IIR filter for this specific problem - finding smoothed first derivative from signal with noisy slope?


texts: Audio Steganography: Inaudible Audio Watermarking for Source Identification

texts:
What is a relatively inaudible audio watermarking technique to overlay a broadband audio signal on top of playing music without analysing the music for low-bitrate source identification? Ideally should be resistant to multi-path, reverb and resonance.

I'm not trying to do any kind of DRM, so I don't care about resampling, time compression or frequency shifting - I just want to figure out in real-time which source is playing without adding too much audible noise to the signal.

At the moment I'm broadcasting a low bitrate ultrasonic signal on top of the music, modulated with GMSK around a carrier frequency of 19500Hz, but the high frequencies attenuate fast and the high frequencies are still audible.

I'm considering broadcasting a super wide-band PN-sequence that sounds like pink noise. I only need to sustain a transfer rate of about 3-4 bytes (24-32 bits) per second. But if I stretch out the symbol length, I'm worried about environmental noise.

I could add a delay if I need to pre-analyse the output buffer before, but I hope that I can introduce a very soft full-spectrum "hiss" or hum that can carry my source identification signal.

I don't need to track the source during periods of silence, so there will always be some music to mask my signal.

Update: I've found some prior audio steganography work, but I don't know what the trade-offs are for real-time modulation (esp. on a low-power embedded device) and if this would at all work without look-ahead on the output buffer.


texts:
I'm trying to do global registration algorithm for panorama like described
here. But I don't understand how they deal with false pairs.

Also if expand the system

u1*I*P1=u1*I => P1=I

u2*I*P2=u2*A21*P1 =>P2=A21*P1

u(i) coefficients are reduced? Or I interpret system of equations wrong?


texts: Can you quarter the processing time for real, symmetric FFTs?

texts:
The DFT of a real signal is Hermite-symmetric, so you can roughly halve the computation time/memory by not bothering to calculate half the values of the spectrum (and complex conjugating the existing values and copying them to the second half if needed).  So the rfft operation takes N samples and outputs N/2 spectrum bins in half the time, for instance.

Signals which are even-symmetrical and real have even-symmetrical and real spectra (and real spectra take half the memory to store as complex spectra), so for symmetrical input, can the calculation be halved again by only using the first half of the signal (N/2) to generate the first half of the spectrum (N/2)?  How?

Is there a way to take the regular FFT of the N/2 half-signal and manipulate the output to produce the N/2 half-spectrum?

(real/odd ⇔ imaginary/odd would work, too, but real/even ⇔ real/even case is simpler to follow.)


texts: How is time or sampling rate accounted for in this model?

texts:
A basic model of coupled strings (eg. piano) is provided here as:

The principle is that it has two identical string simulations each formed by a delay line and LPF. The outputs of these are summed at the "bridge" in the center. The sum is multiplied by $-Hb(z)$ to get the bridge output, which is then subtracted from each of the string loops to simulate damping from energy lost to the bridge.

The bridge equations are:

Where $Rb(s)$ is the impedance of the bridge and $R$ is the impedance of the string.

Assuming I got that right, here's what I don't understand:

If you are taking the summed output of the two string sims and multiplying them by a certain constant $Hb(s)$, then subtracting this from the string loops every sample, how is time or sample rate being accounted for? Ie. If you run it at 96 kHz wouldn't you do the subtraction twice as often as if you run it at 48 kHz doubling the damping? Or otherwise at what interval are you supposed to perform the subtraction?
What does it mean when an expression uses $(s)$ vs. $(z)$? Eg. $Hb(s)$ vs. $Hb(z)$?

Thanks.


texts: What function to apply to a signal to expose its burstiness (i.e., make the burst part more evident than the non-bursty part)

texts:
I am working on some signals one of which I have attached below:

What can I do (i.e. functions that I can apply to this signal) to get something like the signal below:

Note how one can draw a line at 0.25 to differentiate between rownumber <13000(or close), and rownumber>13000

-currently, I have tried checking how far each value of the signal is away from the mean, but that didn't seem to help either (i.e., the signal still did not look anywhere near my desired). By the way, I am "not" a signal processing guru, which means I probably would appreciate explanation intended for someone with little experience in signal processing.

Edit: Based on Dan's answer, I decided to clarify something on the question. I am using this signal for a machine learning related research, and the assumption is that I only have the sparse area, but this figure helps give insight into the nature of the dense part of the signal. Hence, in the above paragraph when I said I tried the deviation from the average, I mean as I received new signals (which looks like the dense part shown in the fig), I tried to find out how far it is from the mean of the sparse part of the signal. So I guess my new revised question will be- if you have the sparse part of the signal, and you want to use its nature, shape or amplitude to differentiate between the dense (which at the initial time is unknown) and the sparse signals, how would you do this?


texts:
I'm trying to work out how this table has been calculated in the attachments under part c. I thought it was $h[0]\star x[n]$ on the 4th row (first row with no label) and then $h[1]\star x[n-1]$ on the row below and then $h[2]\star x[x-2]$ but this doesn't get me the same values.

Would anyone be able to tell me?

Thanks


texts: 2D Convolution as a Doubly Block Circulant Matrix Operating on a Vector

texts:
I was reading Fundamental Image Processing, Chapter 5 (Image Transforms), I encountered the following problem:
Given the arrays $x_1(m,n)$ and $x_2(m,n)$ as follows:

Write their convolution $x_3(m,n) = x_2(m,n)\ast x_1(m,n)$ as a doubly block circulant matrix operating on a vector of size 16 and calculate the result.
 Can anyone please explain what the meaning of problem is and how it should be solved?


texts:
I'm working on the processing of raw FMCW radar data.

In my case the radar works in the [6.014-7.76] GHz band (I have 187 frequency points starting from 6.014 GHz and ended with 7.76GHz) and the data that I collect are organized in a complex-double numbers (I&Q) matrix [Number of frame x 187] per each couple tx-rx antennas.

Unfortunately the radar company doesn't want to provide me with so much information about how the radar works and I am not much of an expert. They said that the collected data are the channel response in the frequency domain.

Now if I plot one received frame I get the results shown in the figure 1 and if I do 1D fft on the detrended raw data I get the figure 2 ( RTI (range time intensity) plot where you can see two humans moving forth and back in front of the radar).

I'm a little bit confused: I thought that the fft processing were to do on the dechirped time domain data to get the RTI plot but I'm doing it on frequency domain and I'm getting result that makes sense, what am I missing here?


texts: Why samples per symbols $\geq 2$ in the GNURadio Constellation Modulator block?

texts:
Let's say I want to emulate the design of a modulator I am building "in real life".  I'd like the components, input data rate, and output bandwidth to be the same as in the real system.  If I use the GNURadio Constellation Modulator block for this, it imposes a lower limit on samples/second.  This results in a different bandwidth for the GNURadio graph, than the real system.
Why does this block impose this restriction?  And, is there a canonical (or preferred) way around it, when emulating real systems?


texts: Solving LASSO (Basis Pursuit Denoising Form) with LARS

texts:
I'm now working on using LARS (Least Angle Regression) algorithm to solve a LASSO problem in Basis Pursuit Denoising form like:

\begin{align*}
\quad && \arg \min_{\beta}{\left\| y - X\beta \right\|}_{2}^{2} + \lambda {\left\| \beta \right\|}_{1} && (1)
\end{align*}

I have done many search jobs on the web, what confused me is LARS is more like an algorithm to solve problems in this form:

\begin{align*}
\quad & \arg \min_{\beta} {\left\| y - X\beta \right\|}_{2}^{2}
&& \text{subject to.} && {\left\| \beta\right\|}_{1} \leq t && (2)
\end{align*}

I know these two forms is equal to each other in theory and LARS can solve problem(1) for all 0<λ≤∞.
But I want know how to use LARS to solve problem(1) in practice, which means how can I write a function like:

function lars(X, y, lambda){
 ...
 return beta;
}


texts: Causal system, order of numerator and denominator

texts:
Is it true that if the order of numerator of a system function is larger than its denominator, then we can say that the system is causal? (Without finding R.O.C.)

Is there similar thing for Stability?


texts:
Folks I've written a small utility that extracts some (and only some) header information from Tifs.

To test it fully I need some Big Endian Tifs, but I cant find any, and I can't find anything to generate them.

Does anyone know of a public repository for tifs with big endian encoding? Or of something that can generate some?

Thanks, apologies if this is off topic.


texts: Weiner filter vs Analytical Linear optimal estimator

texts:
I've been learning recently about estimation of stochastic processes in general, and linear estimation in particular.

On one hand I've seen the following estimation of $X[n]$ from $Y[n]$ -

$X_{opt} = \frac{COV(X[n],Y[n])}{VAR(Y[n])}(Y-E[Y]) + E[X]$

On the other hand we have Weiner filtering, for which -

$X_{opt-w}=h[n] * Y[n] \space\space;\space\space   H(e^{jw})=\frac{S_{XY}(w)}
{S_{YY}(w)} \space\space ; S_{XY},S_{YY}:= PSD \space\space; H(e^{jw})\leftrightarrow h[n]$

All in all, both claim to be optimal linear estimation of $X[n]$ from $Y[n]$, I've been playing with few examples where I could show that $X_{opt-w}=X_{opt}$, but couldn't really understand if that's always the case.

My questions are -
What is a good rule of thumb to prefer each one of them? Will i want to choose $X_{opt}$ only in the case that $E[X[n]],\space E[Y[n]] \ne 0$? When it applies that $X_{opt-w}=X_{opt}$?

Thank you very much!


texts: Impact of the WSOLA analysis increment on the time stretching factor

texts: WSOLA, basically, is the same as OLA but with a dynamic hop size. That part is clear to me. What is not so clear is the following: The stretching factor is determined by alpha = Hs / Ha and WSOLA uses a small increment, obtained through cross-correlation, to determine where the next input frame will be extracted. If that's the case, doesn't that small increment introduce a deviation between the requested alpha and the real one? After all, that increment determines a new analysis hop size, Ha' = Ha + incr, meaning that, if we maintain the same synthesis hop size, the alpha will change.


texts: Nonuniform sampling points based on the first order derivative of signal amplitude

texts:
Problem:

I need to measure a experimental data set $y(x)$.

However, the experiment is extremely time-consuming and I don't have the luxury to have linearly spaced sampling point $x$.

I need the density of sample points to be proportional to the first order derivative of the data $y^{'}(x)$.

Progress:

first measure $y(x)$ once with uniform x as a test run, then compute $y^{'}(x)$.
Then fit $y^{'}(x)$ to $pdf(x,\alpha, \beta, ... )$.
Then I use the parameters to calculate new data points with inverse CDF $X=icdf[linspace[pdf(min),pdf(max),N],\alpha, \beta, ...]$

The problem is, this method only apply to data sets that is distributed in common probability distribution functions.
Is there any way to make this work with arbitrary data ?
Thanks in advance !


texts: Sine Wave spectrum showing multiple frequencies - Audacity

texts:
I am generating a sine wave using audacity
with 2000 Hz, 30s duration, amplitude 0.5

48 kHz is the sampling rate of the recording device.

I just play and record the generated sound, this is the spectrogram that I get. (Top is original while bottom is recorded) Frequencies other than 2kHz are also showing up in the spectrogram with significant power density. (Using 2048 point hanning window for FFT).

This is the time domain representation for the generated and recorded signal, this is shown to ensure clipping is not happening.

Why does the spectrogram show these higher frequencies and how could I remove them.

I could not find answer to my query on other similar posts in DSP stack exchange thus am creating this post.


texts: This is the signal whose FT i need to find, at first i thought that i could solve this as a convolution of two rectangular pulses, but i could not find pulses that fit into this (it turns out that parts from -T/2 to 0 and T to 3T/2 are always somehow narrow). I could also try by definition but that's just too much work and very time consuming, the third option is derivative theorem, but i don't know how to use it, since i never actually used it. Any help appreciated!


texts: Why do the DTFT and FFT give me completely different results for magnitude at a specific frequency?

texts:
I am trying to write a program to compute the magnitude and phase of a specific, non-integer frequency component (i.e. given a sampled finite signal of length $N$, I want to know the magnitude and phase of the spectrum at frequency 40.678 Hz).
I tried using the equation for the Discrete Time Fourier Transform:
$$
X(\omega) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n} = \sum_{n=0}^{N-1} x[n] e^{-j\omega n}
$$
such that
$$
X(2\pi\cdot40.678) = \sum_{n=0}^{N-1} x[n] e^{-j 2\pi\cdot40.678\cdot n}
$$
I did this as follows:
t = np.linspace(0,1,100)
x = np.sin(2*np.pi*40.678*t)

Xf = 0
for n in range(len(x)):
    Xf += (x[n]*np.exp((-np.complex(0,1))*2*np.pi*40.678*n))
print(abs(Xf))

Which printed $0.6149015687259044$.
I know that this is wrong because if I do an FFT as follows:
X = np.fft.fft(x)
plt.figure();
plt.plot(abs(X))

I can see peaks of magnitude $50$ at around (what I assume after correcting the axis is) 40.678 Hz.

What is the reason for this difference? I also tried the equation of the DFT just like I did with the DTFT code above this time to measure magnitude at 50 Hz of a sinusoid of frequency 50 Hz and I got a completely wrong value again as compared to what the plot of the FFT shows.
What causes this discrepancy?
Thank you.


texts:
I came across this process when working on OFDM in wireless communication. It's mentioned if we have two copies of a signal resulted from SIMO system (Single Rx Multiple Tx), we can process the noise whitening, after the FFT, in order to improve the quality of the signal.

What does mean the noise whitening, and what's its process?

Here is an example,  in the paper.

After taking the DFT in OFDM, we supposed to get frequency response, so how w[k] is obtained?


texts: Deciding on the FFT length for water waves in flume (wave facility)

texts:
Soon I'll be conducting experiments in a wave facility (flume) for my MSc thesis. And I would like to analyse wave spectra (amplitude and energy density).

The time step (or 'accuracy') of the wave maker, is by default dt = 0.001 s.

The sampling frequency of the measurement instruments (wave gauges) is 100 Hz (previously 25 Hz).

I'll start with bichromatic waves, which can be more or less be treated as being deterministic. I would like to prevent spectral leakage as much as possible.

The fundamental frequency is different for different experiments: short waves with a wave period of around 2.5 s. The wave maker includes bound waves (subharmonics), with the wave period Tbound one over the difference frequency of the short waves in the bichromatic wave train.

So far, previous students have set the FFT length to a power of 2. From what I read this is not necessary with modern implementations of the FFT.

One example of a previous student:

Fs = 25 Hz,

N = D*Fs = 2^13 (power of 2)

This results in a duration of the experiment of
D = N/Fs = 327.68 s

with a spectral resolution of
df = 1/D = 0.00305175781 Hz.

To avoid spectral leakage, frequencies of all my wave components should be a multiple of this number if I am correct.

Since I've set my wave gauges at a sampling frequency 100 Hz, I'm wondering if it's fine to increase my duration to D = 1000 s, so N = 100*1000 and is not a power of 2.

The resulting spectral resolution is:
df = 1/D = 0.001 Hz
and I'll round the frequencies of all my wave components to a multiple of this number.


texts: Design of the low pass filter when decimating

texts:
I am filtering an audio signal into various sub-bands in C and would like to decimate the signal by 2 in the lower frequency bands to reduce processing overhead. I get that this involves low-pass filtering the signal and then throwing away every other sample, but how do I go about designing the low-pass filter optimally?

At the moment I have a signal sampled at 44.1 kHz. Say that for signals in bands limited below 8 kHz, I low-pass filter with a 2nd order butterworth filter with a cutoff of 9 kHz and then decimate. This works fine as I filter well below the new Nyquist limit of 44100/4, but it feels very unscientific.

What is the best filter type for this purpose. Looking at the matlab documentation, the decimate function uses an 8th order Cheby Type 1 filter by default, but I have no idea why. Is there some rule-of-thumb that I'm missing?


texts: wavelet analysis normalized to have unit energy

texts:
I am analyzing some time series of sea surface water temperature and am using wavelet analysis to do so. I am computing the continuous wavelet transform and then removing some specific frequencies and then converting the series back from the frequency to the time domain. I am following the methods of torrence and compo 1998 (Practical guide to wavelet analysis) to perform this operation but am a bit confused as to why my results are not as expected. As an example, I will fist decompose the time series into the frequency domain and then reconstruct the signal without altering any of the coefficients:

In the attached figure we have the original series in red and the reconstructed series in blue, straight away we can see that the the reconstructed signal is much lower than the original series. The main reason I think this occurs is that the wavelet transform normalizes the series to have unit energy (in order to compare the coefficients across scales). Therefore, I was wondering if it was possible to convert the reconstructed signal back to have the same magnitude (?) as the original series i.e. reverse the normalizing to have unit energy stage?


texts: Is non-causal, non-LTI solution to difference equation correct?

texts:
Find a solution to the difference equation $y[n]-\frac{5}{6}y[n-1]+\frac{1}{6}y[n-2]=\frac{1}{3}x[n-1]$ that is neither casual nor LTI, where $y[0]=y[1]=1$ and $x[n]=\delta[n]$

The homogenous solution for the difference equation is

$$y_h[n] =A_1(\frac{1}{2})^n + A_2(\frac{1}{3})^n$$

Solving for specific solution ($y_h[0]=y_h[1]=1$):

\begin{equation}
y_1[n]=4(\frac{1}{2})^n-3(\frac{1}{3})^n \tag{1}
\end{equation}

We want to "modify" $y_1[1]$ or $y_1[0]$ or $y_1[-1]$ because the right hand side of the difference equation is non-zero at $n=1$.

Since $y[1]$ and $y[0]$ are specifically set to $1$ in the question, we must change $y[-1]$.

$$y[1]-\frac{5}{6}y[0]+\frac{1}{6}y[-1]=\frac{1}{3}x[0]$$
$$1-\frac{5}{6}+\frac{1}{6}y[-1]=\frac{1}{3}$$
$$y[-1]=1$$

Since $y_1[-1]=-1$, $y_1$ will be the correct output only for $n \geq 0$. This means that we need to find another "piece" $y_2$ for $n < 0$ such that $y_2[0]=y_1[0]$.

Solving another specific solution with $y_h[-1]=y_h[0]=1$:

\begin{equation}
y_2[n]=2(\frac{1}{2})^n-(\frac{1}{3})^n \tag 2
\end{equation}

Hence, our solution is:

$$
y[n] =
\left\{
\begin{array}{ll}
       4(\frac{1}{2})^n-3(\frac{1}{3})^n & n\geq 0 \\
      2(\frac{1}{2})^n-(\frac{1}{3})^n & n \lt 0 \\
\end{array}
\right.
$$

which can be written as

$$y[n]=4(\frac{1}{2})^n-3(\frac{1}{3})^n-2(\frac{1}{2})^nu[-n-1]+2(\frac{1}{3})^nu[-n-1]$$

Is there a faster way to approach this problem (possibly using Fourier transform)? How would you solve this for complicated inputs, e.g. $x[n]=u[n]$?


texts: I'm trying to develop an automated state-detection algorithm, basically expanding on MatLab's "statelevels" function to multiple (i.e. 3+) states, however I've hit a brick wall in trying to figure out how to define histogram bin edges without prior knowledge of the signal and have not been able to find any literature which isn't locked behind a paywall. Does anyone have suggestions, either for an algorithm to achieve multi-state detection or for publicly-available literature?


texts: How do I get the Time information after STFT?

texts:
I'm tring to code a STFT program in Matlab. I got results but I'm confused a bit.
I have a audio signal of 900000 pts sampled at 1000Hz (15 min signal).
I've taken fft with following parameters =>

hop size of H = 1400 pts
the window length M = H * 6

So, 5/6 of the window overlaps with the next one.
The program runs for 637 iterations, as (637*1400 + M = 900200)

Now, I get a 637 pt vector with frequency in them. I should get a frequency vs time plot. But it only has 637 frequency bins, for a signal of 900000 time sample.

How, can I find the change of frequency with Time (900000 sample) or
15 min?
Moreover, how does changing the window length and overlap affect the outcome of STFT? The code is below.

%hop size
L = 1400;

%window length & Formation
D = 6;                        %User Defined Parameter in Cooper's paper
M = L*D;
w = hanning(M, 'periodic');    %window formation

k = 1; % k =index
j = 1; % j =result vector index
while k + M <= xlen

  %data chunk
  %x = x(k:1400);

  %window length & Window
  xw = x(k+1:k+M) .*w;                  %windowed signal

  %zero padding
  b = 5;                           %Zero Padding Factor 'b'
  xw = [xw; zeros(M*5,1)];
  [nft,c] = size(xw);

  %plot(abs(fft(xw)))
  X = abs(fft(xw));
  XF = X(1:length(X)/2);
  [q,ind] = max(XF);
  v(j) = ((fs/2)/(nft/2)) * (ind-1);  %convert Bin to frequency
  k = k+L;
  j = j+1;
end

figure(1);
plot(v)
axis([0 1400 0 100])

figure(2);   %zoom version
plot(v)


texts: In OFDM, what is the length of the cyclic suffix that is added before windowing?

texts: I searched for length of cyclic suffix but all I can find information about is the cyclic prefix, I want to know also if there is short guard interval and long guard interval does it affect the length of the suffix


texts: Can Log Gabor Filters be used on an edge detected image

texts: i'm doing a project on log gabor filters, can i use log gabor filter on an edge detected image for feature extraction,and use these features for classification by ANN


texts: Non-symmetrical DFT of real valued signal

texts:
I have a homework assignment to take the $64$ DFT of $\cos((5\pi/32)n)$ and $\cos((5\pi/64)n)$. However, for the second case I obtained DFT spectrum such that I have that famous symmetry in real part and imaginary part of it because it is real valued signal. However, for the first signal I found something that has no symmetry in imaginary part of the DFT, but it's real too.

Why there is no symmetry although signal is real? Does it have anything to do with "windowing" the signal?

I wanted to calculate it analytically so that here is my code and the results:

....
t = 0:63;
x1 = cos(5*(pi/32)*t);
X1_analytic = zeros(1,64);
temp = 0;
for k = 0:63
    for n = 0:63
        temp = temp + x1(n+1)*exp(-1i*2*pi/64*k*n);
    end
    X1_analytic(k+1) = temp;
    temp = 0;
end
subplot(5,1,1);
stem(0:63,abs(X1_analytic),'--*r');
subplot(5,1,2);
stem(0:63,(X1_analytic+conj(X1_analytic))*(1/2));
title('Real part');
subplot(5,1,3);
stem(0:63,imag((X1_analytic)));
title('Imaginary part');
subplot(5,1,4);
stem(0:63,abs(X1_analytic));
title('Abs value ');
subplot(5,1,5);
stem(0:63,angle(X1_analytic));
title('Phase ');


texts: Periodicity of Fourier series's coefficients

texts: Why exactly continious Fourier-series's coefficients aren't periodic like coefficients of discrete Fourier series (DFS)? $e^{-j2\pi}$ is periodic in both sequences.


texts:
I like to work with a convolutional neural network in combination with a spectrogram as input.

Assuming the spectrogram has the dimensions $T\times F$ (time$\times$frequency). Is it more natural to see the spectrogram as a 1D image with $F$ channels and the size $T\times 1$ or as a 2D image with the size $T\times F$?

Depending on this I work with 1D convolutions or with 2D convolutions. But I'm not sure whats the more natural way.

Thank you very much


texts:
Assume that I have L points and want to calculate its Fourier transform with N points (N is much larger than L, link N=256, while L=10), I am think if there is way to calculate it with lowest complexity.

At first, I try to calculate it by the definition Fourier transform, according to the definition, the complexity cost should be 2NL.

Then I try to calculate with fast Fourier transform (FFT). As far as I understand, I should make N points input (L points with N-L zeros) to produce N points output. In this case, the complexity should still be the NlogN+N, it is small than 2NL.

Thinking that L is much smaller that N, the result is unaccepted that it takes the same computational cost with N points input. I wonder do you guys have better idea?


texts: Real time numerical differentiation of signals

texts:
I am trying to implement a PID controller in Python and I am having some problems with real time numerical differentiation of my my discrete signal.

I am using the following method:

$$d{\frac{x_n}{t}} = \frac{x_{n-1}-x_n}{d_t}$$

Where $d_t$ is the time difference of the measurement time of $x_n$ and $x_{n-1}$.

The differentiation is accurate to some point but the results on real time systems are not as accurate as I would like to be. The differential controller does not really improve the stability of the system in most cases. At high gains it starts introducing noise into the system.

Most numerical differentiation methods recommend going into the future ($x _{n+1}$) and sadly for my system that's not possible. Such as:

$$\frac{x_{n+1}-x_{n-1}}{2\times d_t}$$

Are there any methods for calculating real-time differentiation with better accuracy?


texts:
I am still trying to iron out some ambiguities in my understanding of the IFFT component of OFDM modulation schemes.

So here we have a QAM symbol $s_0$ being multiplied with the subcarrier for that frequency bin, $s_1$ being multiplied with carrier 1, $s_2$ carrier 2, and it shows clearly that the sum of these products produces the OFDM symbol in the time domain. This would logically follow the equation:

$
c\left(t\right) = \sum^{N-1}_{n=0}{x_n e^{-i\frac{2\pi n}{N}}}
$

Where c(t) is the output symbol, x is the list of complex QAM symbols and n is the subcarrier index, N is the number of subcarriers. (Note there is no 1/N at the start and i'm not sure why that needs to be introduced)

Why is it that IDFT is stated on Wikipedia to be

$
x_n = \frac{1}{N} \sum^{N-1}_{k=0}{X_k . e^{i 2 \pi k n / N}}
$

Where x is the output list and xn is merely a single output sampling point in the time domain.

Here is an example of application of the DFT/IDFT formula :

Let :
$N = 4$
and :
$
x = \left(\begin{array}{c}x_0\\ x_1\\ x_2\\ x_3\\
\end{array}
\right)
=
\left(\begin{array}{c}1\\ 2-i\\ -i\\ -1+2i\\
\end{array}
\right)
$

then

$X_0 = e^{-i 2 \pi 0.0/4} + e^{-i 2 \pi 0.1/4}.\left(2-i\right) + e^{-i 2 \pi 0.2/4}.\left(-i\right)+ e^{-i 2 \pi 0.3/4}.\left(-1+2i\right) = 2$
$X_1 = e^{-i 2 \pi 1.0/4} + e^{-i 2 \pi 1.1/4}.\left(2-i\right) + e^{-i 2 \pi 1.2/4}.\left(-i\right)+ e^{-i 2 \pi 1.3/4}.\left(-1+2i\right) = -2-2i$
$X_2 = e^{-i 2 \pi 2.0/4} + e^{-i 2 \pi 2.1/4}.\left(2-i\right) + e^{-i 2 \pi 2.2/4}.\left(-i\right)+ e^{-i 2 \pi 2.3/4}.\left(-1+2i\right) = -2i$
$X_3 = e^{-i 2 \pi 3.0/4} + e^{-i 2 \pi 3.1/4}.\left(2-i\right) + e^{-i 2 \pi 3.2/4}.\left(-i\right)+ e^{-i 2 \pi 3.3/4}.\left(-1+2i\right) = 4+4i$

$
X =
\left(\begin{array}{c}X_0\\ X_1\\ X_2\\ X_3\\
\end{array}
\right)
=
\left(\begin{array}{c}2\\ -2-2i\\ -2i\\ 4+4i\\
\end{array}
\right)
$

If we pretend this worked example is IDFT and not DFT, it clearly shows that for each output sampling point, Xk, in the time domain, every point in the constellation was multiplied by a single frequency bin and then summated. This is directly backed up by the fact that IDFT is stated to have a complexity of O(n^2) (as far as I'm aware IFDT and IFFT are the same equation but implemented differently algorithmically).

Why is this? This directly contradicts many OFDM examples and diagrams I've seen and it fundamentally makes sense to me that a single frequency would be phase and amplitude shifted by a single complex input from the constellation, and not by all of them.


texts:
I have this image

I would like to remove the shadow from the image. I know a lot of different methods like certain morphological operations have been used to remove shadows:

I have created this mask for the same image

Are there some other methods I could try using this mask that I have created?

EDIT:

input image and mask of same size as requested :

EDIT 2: i generated a 1D invariant image however its not perfect

  I = imread('shadow.jpg');
       J = im2double(I);

      R = J(:,:,1);
      G = J(:,:,2);
      B = J(:,:,3);

     [len,wid] = size(R);

     % Generation of 2-D Log Chromaticity Image.
     for i = 1:len
        for j = 1:wid
           if ((R(i,j)*G(i,j)*B(i,j))~= 0)
              c1(i,j) = R(i,j)/((R(i,j)*G(i,j)*B(i,j))^(1/3));
              c2(i,j) = G(i,j)/((R(i,j)*G(i,j)*B(i,j))^(1/3));
              c3(i,j) = B(i,j)/((R(i,j)*G(i,j)*B(i,j))^(1/3));
           else
              c1(i,j) = 1;
              c2(i,j) = 1;
              c3(i,j) = 1;
        end
    end
end

rho1 = mat2gray(log(c1));
rho2 = mat2gray(log(c2));
rho3 = mat2gray(log(c3));

X1 = mat2gray(rho1*1/(sqrt(2)) - rho2*1/(sqrt(2)));                                         %(1/sqrt(2); -1/sqrt(2); 0)
X2 = mat2gray(rho1*1/(sqrt(6)) + rho2*1/(sqrt(6)) - rho3*2/(sqrt(6)));   %(1/sqrt(6); 1/sqrt(6); -2/sqrt(6))

theta = 120;

InvariantImage = cos(theta*pi/180)*X1 + sin(theta*pi/180)*X2;
imagesc(InvariantImage); colormap(gray)

Cannot understand what i am doing wrong out here please help?


texts: State space equation from differential equation

texts:
I have very general system. I don't know whether it is electrical or mechanical or whatever. This system can be modeled by the following differential equation

$$\dot q = \frac{Tf_1-f_2}{T+1}$$

where:

$\dot q$ would be equivalent to current in an electrical system
$f_1$ and $f_2$ would be equivalent to current sources
$T$ is a given constant

From what I see, if $f_1$ and $f_2$ are constant, there is no way for $\dot q$ to change. How is it possible that this system has a state space representation? And how do I get this state space representation?


texts: What is the meaning of the DTFT of the unit impulse sequence?

texts:
In an exercice, I'm asked to draw the $X_{imp}(\omega)$ Discrete-Time Fourier Transform (DTFT) of the $x_{imp}(n)$ unit impulse sequence defined as:
$$
x_{imp}(n) = \begin{cases}
  1 & \text{if } n = 0 \\
  0 & \text{otherwise}
\end{cases}
$$

Thanks to another question, I think I understand the difference between the DFT and the DTFT.

So, given $X(\omega) = \sum_{-\infty}^{\infty} x(n)\cdot{\rm e}^{-j\omega n}$, I end up with $X_{imp}(\omega) = 1$ because all terms of the summation would be $0$, except at $n=0$ where $x(n)\cdot{\rm e}^{-j\omega n} = 1$.

If my answer is correct, does that mean the unit impulse sequence would have the same amount of "energy" (← is this the right word?) for all frequencies, including audio, radio frequencies X-Rays, light or whatever else. That seems somewhat irrealistic to me. So, where is the flaw in my reasoning?


texts:
In ofdm systems after applying IDFT we take the last N/x samples and place them before these N samples. This is called cyclic prefix(CP).So why to use a cyclic prefix?

Why is the CP known to combat the  inter symbol interference.Why can't we use some junk data or zeros over there? Can't this junk data combat the ISI ?.

Regards,
Phani Tej


texts: Question about ramp filter used in filtered backprojection

texts:
Question is this.
First, a ramp filter (in frequency domain) is defined by $H(Q)=|Q|$.
What are the responses of a ramp filter to
(1) a constant function $f(r)=c$ and (2) a sinusoid function $f(r)=\sin(wr)$?
What does the response mean? Following is my work.

My work:

First, take fourier transform of a function $f(r)=c$. It is
$\int_{-\infty}^{\infty}f(r)e^{-2i\pi rQ}dr=c\delta(Q)$. Then multiply ramp filter and take inverse fourier transform. It is $\int_{-\infty}^{\infty}c\delta(Q)|Q|e^{2i\pi irQ}dQ=0$??
Similarly, $\int_{-\infty}^{\infty}\sin(wr)e^{-2i\pi rQ}dr=\frac{\delta(Q-w/2\pi)-\delta(Q+w/2\pi)}{2i}$. So Applying the ramp filter and i.f.t gives
$\frac{w(e^{iwr}-e^{-iwr})}{4i\pi}=\frac{w\sin(wr)}{2\pi}$.

It this right?


texts:
I am doing the spectral analysis of some antenna recorded signals by computing PSD (logarithmic absolute value of the FFT) and the spectrogram of the signal. The problem is that I get some unexpected spikes in the signal that appear through all the signal, as it can be seen in the next figure:

The two peaks that appear are the spikes I am speaking about, and they are the same in all the signal. The center spike is not located in $0\textrm{ Hz}$, so it is not a DC offset of the signal. This is giving me some problems, so I would like to know what can they be or how can they be removed (maybe I am not computing the spectral analysis correctly).

I have also performed Welch PSD estimations, and the peaks still appear. I have been told that it can be some internal interference form the electronics recording the data or some issue with ADCs.


texts: Impulse response of an infinite-length system

texts:
I have a system which is expressed as following:
$$y[n] = x[n] - \alpha y[n-N]$$ where $y[n]$ is the output and $x[n]$ is the input.
When I take the z-transform of both sides, I've found the transfer function as $$H(z)=\frac{1}{1+\alpha z^{-N}}$$.
However, I need to find the impulse response of this system but I couldn't find $h[n]$. Can anybody help me how to solve this problem?


texts: Digital Communications block diagram question

texts:
I am new to digital communications and trying to understand the basic underlying concepts of digital comm systems. So far , studying the course , I have some questions which if anyone can answer would be great.

In the transmit section, to represent the data (be whatever modulation we are using or No modulation), we have to represent bits/symbols (bits-no modulation, symbols-with modulation) in terms of pulse shaping function.   Am I right?

And mostly we use root raised cosine signal to represent the symbols in order to restrict the leakage and have square(approx) shaped spectrum


texts: Long term spectral flatness measure voice activity detection

texts:
This is the Long term spectral flatness measure code that I have written so far. Can someone please help me, I don't know if I'm writing it correctly. It is  working fine and gives high peak for speech region and sometimes output high flatness for regions where it should actually give high peaks.

nfft = 1024;
s_bin = ceil(nfft*(500/fs))+1;
e_bin = ceil(nfft*(4000/fs));

[S]=spectrogram(signal_zpad,hann(frame_length),frame_shift_length,nfft,fs);

M = 10;
R = 30;
spec_size = size(S,2);

SS(1:spec_size) = 0;
AM(1:e_bin) = 0;
GM(1:e_bin) = 1;
RR(1:spec_size) = 0;

for m=1:spec_size
    for wk=s_bin:e_bin
        for n=m-R+1:m
             if n>0
                c = 0;
                for p=n-M+1:n
                    if p>0
                    SS(n) = SS(n) + abs(S(wk,p))^2; c=c+1;
                    end
                end
                if c<M
                    AM(wk) = AM(wk)+ SS(n)/c;
                    GM(wk) = GM(wk)*(SS(n)/c);
                else
                    AM(wk) = AM(wk)+ SS(n)/M;
                    GM(wk) = GM(wk)*(SS(n)/M);
                end
                SS(n) = 0;
             end
        end
        if m<=R && c>0
            AM(wk) = AM(wk)/m;
            GM(wk) = nthroot(GM(wk),m);
        elseif m>R
            AM(wk) = AM(wk)/R;
            GM(wk) = nthroot(GM(wk),R);
        else
            AM(wk) = 0;GM(wk)=0;
        end
        RR(m) = RR(m) + log10(GM(wk)/AM(wk));

    end

        L(m) = RR(m);

    AM(1:e_bin) = 0;
    GM(1:e_bin) = 1;
    RR(m) = 0;
end


texts:
Generated a 1kHz square wave in Audacity and normalized it to 0dB (fs = 44100Hz):

Added a 24dB/oct high pass filter (std Audacity plug-in), cutoff frequency set to 100Hz. Result:

Zoomed version:

Why "clipping" so badly?
EDIT: When normalized to -3dB before filtering there's no "clipping" present.
Filtering same signal with low pass filter does not have this issue.


texts: What is the normalization factor for the quincunx lattice wavelet transform, and how do you find it?

texts:
On pages 57-60 (preview was available last I checked, images here in case), there is a quincunx lattice transform described.

Lattice:

o • o • o • o •
• o • o • o • o
o • o • o • o •
• o • o • o • o
o • o • o • o •
• o • o • o • o
o • o • o • o •
• o • o • o • o

Basically you do these Predict operations on the black points:

x[ m][n ] -= 1/4 * ( LEFT + RIGHT + DOWN + UP )

Where $LEFT = x[ m ][ n-1 ]$, $RIGHT= x[m][n+1]$, $DOWN=x[m+1][n]$, $UP=x[m-1][n]$.

You then do updates on the white points:

x[ m][n] += 1/8 * ( LEFT + RIGHT + DOWN + UP )

Then you will never touch the black values again, so you effectively have:

o x o x o x o x
x o x o x o x o
o x o x o x o x
x o x o x o x o
o x o x o x o x
x o x o x o x o
o x o x o x o x
x o x o x o x o

You turn your head 45 degrees to see this is just another rectangular lattice, and you label them odd/even again:

o   o   o   o
  •   •   •   •
o   o   o   o
  •   •   •   •
o   o   o   o
  •   •   •   •
o   o   o   o
  •   •   •   •

You repeat this again and again, until you have 1 "average" left.

Now in the Haar wavelet transform, there is a power loss in each level that we correct with a normalization factor of √2.

Here, there is a computed power loss factor of about 1.4629 after the first step of the first level (found by running 5,000,000 transforms on random data and finding the ratio of powerBefore/powerAfter and averaging).

I don't know how to show/compute how this power loss is found, and where the 1.46 number comes from.


texts: What's the rationale behind this pulse shaping function?

texts:
Section 13.2 of the HD Radio standard for AM defines a pulse shaping function that appears to be some strange variant on a raised cosine filter:

More specifically, H(ξ) is a raised cosine filter with total width (1+α)T, and W(ξ) appears to be a "stretched" version. Why do they do this?


texts: How to reliably compute the group delay of a comb filter

texts:
I applied the following FIR comb filter in real-time:

y[n]=x[n]-x[n-40]

Since this is an FIR, the group delay is D=(N-1)/2=20 samples. After applying the filter to a signal, I tried to use cross correlation between the filtered and unfiltered signal, to reproduce D computationally by determining the argmax of the cross correlation (I do have a need to the delay this way). The issue is that I get too peaks in the cross correlation, one at zero lag and another at 20 lag. But the peak at zero lag is the maxima which means the peak at 20 lag which is the correct lag is ignored. This method work really well with other filters like averaging filters.

Does anyone know while I get a the peak at zero which is overshadowing the real peak? Is this normal for comb filters? Is there another method to compute delays using the filtered and unfiltered signal other than cross correlation?


texts: CDMA scheme that preserves constant envelope (NRZ/ASK shape)

texts:
(This is a followup question to Is this a "valid" CDMA modulation?)

Is it possible to find a CDMA scheme in which the receiver can only have a simple slicer?
In other words, for the receiver, the signal should look like an ASK-like signal, i.e. it can detect a bit of 1 (above the average) or a bit of 0 (below the average) but cannot detect signal levels. It also cannot perform despreading in analog domain.

Further explanations:

It seems to me that any CDMA scheme that uses Walsh codes (or PN codes) result inevitably in a signal with multiple levels when the transmissions of multiple stations are summed. For example, even in the simple example from https://en.wikipedia.org/wiki/Code-division_multiple_access#Example, the output signal is $(0, −2, −2, 0, 2, 0, 2, 0)$, consisting of three different signal levels. For demodulation, the receiver needs to sample the signal levels in order to sum them together.

On the contrary, operating in $GF(2)$ seems to preserve the desired signal properties. However, it seems real spectrum spreading is not possible with this method, since only one station would transmit a bit at a given time (and while technically being CDMA, being more time division multiplex, see Is this a "valid" CDMA modulation?)


texts: How to crop a song and mantain the same values in the sample vector?

texts:
Does anyone know how can I crop an audio signal and maintain the same values in the sample vector?

I am working on a audio signal processing project for song recognition in Python and I wanted to crop a song in Audacity. To my surprise, the audio samples (for both the left and right channel) in the cropped version were "corrupted", meaning that they were different from the corresponding samples from the original song.


texts: How do we impose Hermitian Symmetry in optical OFDM?

texts:
I’m currently trying to wrap my head around optical OFDM. I’m reading various tutorials and articles, but I’m hitting a stumbling block trying to understand Hermitian Symmetry and how it is applied.
I’m taking the block diagram of an OFDM system from Figure 3 in this paper, and working through it step by step starting from the input to the transmitter. It’s probably best that I go step-by-step up until the point where I get confused.

Step 1 - serial to parallel
For the sake of simplicity I’m skipping the “Coding” and “Interleaving” steps from that paper. Let’s start with a serial bit stream we want to send: 1, 1, -1, -1, 1, 1, 1, -1, 1, -1, -1, -1, -1, 1, -1, -1, -1, 1…etc.
We will map this to parallel bit streams, one for each subcarrier. I’m using 4 subcarriers in this example, C1 to C4. The parallel streams are now:

C1
C2
C3
C4

1
1
-1
-1

1
1
1
-1

1
-1
-1
-1

-1
1
-1
-1

-1
1
1
-1

-1
-1
1
1

Step 2 - QAM modulation
Modulate each subcarrier. Here I'll use 4-QAM. Taking 2 bits at a time from each of the subcarrier columns I get:

C1
C2
C3
C4

1 +1j
1 +1j
-1 +1j
-1 -1j

1 -1j
-1 +1j
-1 -1j
-1-1j

-1 -1j
1 -1j
1 +1j
-1 +1j

Step 3 - IFFT
This is where my understanding starts to break down. I understand for “normal” OFDM we’d pass these values to the IFFT to generate the time-domain OFDM signal. However, what exactly are we sending to the IFFT block? Do we send it each row of that table in turn, i.e. the values of each subcarrier at a particular sample time?
Step 4 - Hermitian Symmetry
This is where my understanding really breaks down! I understand that Optical-OFDM require purely real-valued time-domain OFDM signals. I also follow that if the input to the IFFT block has the property of "Hermitian symmetry", then the IFFT output will be the desired real signal.
Hermitian symmetry is "imposed" on Xk symbols by ensuring that:
X0 = XN/2 = 0 and
Xk = X*N-k for 0 < k < N/2
Where N is the IFFT length and k is the carrier number (I assume?)
My confusion here is if we have to "impose" that certain carriers must take a given value (e.g. carrier 0  and carrier N/2 must equal 0) then isn't that somehow "overwriting" the value for that carrier that we have applied at the QAM modulation stage? That seems like it must be wrong to me, but I can't understand why!
So, apologies for the long-winded question, but just to clarify I'm asking if my understanding of the IFFT at Step 3 is correct, and what exactly we are doing when we "impose" Hermitian Symmetry onto our signal prior to the IFFT in Step 4.


texts: How to account for number of electrodes in EEG data while converting to frequency domain?

texts:
I have raw data from 64 electrode EEG in CSV format where there are multiple rows and 64 columns. I am using the following Python code to separate the data into frequency bands:

import numpy as np
import pandas as pd
import csv

# Sampling rate = 500 Hz
fs = 500

# Define EEG bands
eeg_bands = {'Delta':  (0, 4),
             'Theta':  (4, 8),
             'Alpha':  (8, 14),
             'Beta':   (14, 30),
             'Gamma':  (30, 45)}

f = open('output-stddev.csv', 'w', newline='')
w = csv.DictWriter(f, eeg_bands.keys())

data = np.genfromtxt('rawCSV.csv', delimiter=',')

# Get real amplitudes of FFT (only in postive frequencies)
fft_vals = np.absolute(np.fft.rfft(data))

# Get frequencies for amplitudes in Hz
fft_freq = np.fft.rfftfreq(len(data), 1.0/fs)

# Take the mean of the fft amplitude for each EEG band
eeg_band_fft = dict()
for band in eeg_bands:
    freq_ix = np.where((fft_freq >= eeg_bands[band][0]) &
                       (fft_freq < eeg_bands[band][1]))[0]
    eeg_band_fft[band] = np.std(fft_vals[freq_ix])

w.writerow(eeg_band_fft)
f.close()

I believe the code is not correctly thought-out because it doesn't account for the fact that there are 64 electrodes anywhere in it. How do I accurately convert raw EEG data from all the electrodes to the frequency domain?


texts: Explaining / Communicating spectra: Generating Visualizations

texts:
I'm very often in the situation that I need to explain / visualize spectrum of signals or frequency behaviour of systems (esp. filters). My favourite, very "textbook" styled way of communicating these would be, in simpler cases:

¹

Typically, I'd have to visualize something like aliasing, effects of suppressing carriers etc. This also leads to less simple situations, like

¹

Now, drawing these kinds of diagrams is really a burden if done with my vector program of choice (that being Inkscape) – I want to parametrically specify the spectrum, not manually shift, copy, mirror and repeat lines!

So:

How do you approach theses kind of visualizations? Vector graphics programs & patience? Pen & paper? Or do you have a private set of scripts that e.g. draw the diagram with the help of matplotlib or tikz?

¹ Bellanger: Digital Processing of Signals: Theory and Practice. 2nd Edition, Wiley 1988. P. 17 & p. 18.


texts:
In order to submit a document to an online website, they require the following:

Scan the US letter-size hard-copy document at 600 dpi in greyscale
Provide a JPEG, JPG, or PNG file of the scan with a maximum size of 5 MB

When I scan at 600dpi and create a JPG with no compression, I get an 8.5 MB file.  So, I used (FastStone image) software to compress it down below 5 MB.  When I submitted that document, the website automatically responded "document unacceptable - less than 600 dpi."
Is it possible to create what they are requesting?  Any advice would be appreciated.


texts: Synchronizing signals with different sample rates

texts:
I am relatively new to signal processing, so I am still missing some basic knowledge, but I am trying to improve as much as I can. I am trying to synchronize two signals with matlab using cross-correlation and I followed the matlab documentation. The problem is that my two signals are with different sample rates. The first signal is sampled at 1500 Hz and the second one at 800 Hz. Using matlab I downsampled the first signal to 800 Hz - signal1 = resample(signal1, 800, 1500); and the synchronization worked.

Is there a way to avoid downsampling the first signal, because from what I know, when I downsample a signal, I lose information and I would like to avoid that. Furthermore, I am still not sure how to check if it meets the nyquist criteria, which if I have understood correctly, the nyquist frequency must be twice the sample rate. I did a fft on my signal and if I am not wrong the frequency is roughly 750 Hz, so downsampling the first signal is not my best option. I will include a photo of the fft results.

The fft was performed by following the documentation from matlab on how to plot the Single-Sided Amplitude Spectrum of X(t).


texts: How can a CNN account for spectro-temporal constraints in neural data?

texts:
What are there the best ways to leverage the unique "geometrical" constraints of spectro-temporal signal representations (architecture, filter shapes, data augmentation, etc.)? For example, a learned feature in a spectrogram should intuitively be horizontally (temporally) invariant since the window of time when it was captured is arbitrary. In contrast, the specific vertical position of a feature is extremely important since the precise frequency at which it occurred is significant. My research group is trying do design a convolutional neural network (using keras/tensorflow) that takes spectrograms representing neural data and predicts a particular disease biomarker. All of our trained models so far severely overfit and we have somewhat limited data.
Just for further background:
We decided to use CNNs because the visual representations of spectral data are known to be useful for our task.
Even though our data set is not huge, we were hoping that CNNs could be more interpretable in terms of exactly what spectral features are being learned than other ML methods like SVMs or LSTMs.
The models have been regressors so far, but we are open to switching to classifiers if they perform better in the end.


texts: Replicate MATLAB's `conv2()` in Frequency Domain

texts:
When conv2d is on valid mode, the image needs no padding, because the result is the same size as the image.
When conv2d is on full mode, the result is (image_width + kernel_width -1) * (image_height + kernel_height -1).
Then how do I pad the image?


texts: Deep Learning based NDA Channel Estimation

texts:
I was wondering if it is possible to use Deep Learning to estimate the Channel Impulse Response for NDA synchronization?

I understand Deep Learning is not normally used for regression problems but I believe you can approximate it as a classification problem if you make your channel tap magnitudes discrete.

I'm not completely sure how NN would be built since you would have the number of SoftMax functions equal to the number of channel taps along with each channel tap having a set of classification values it could take on.

To remove phase and frequency offset dependence you could take the magnitude of the received complex signal and as an assumption assume the tap with the greatest magnitude is real and positive.


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

texts:
The Short Time Fourier Transform (STFT) is used to identify time localized frequency content of a signal. The STFT operates by chunking an input signal into blocks and performing FFT on the block, and then recording the frequency coefficients for that point in time. However, the STFT enforces fixed time and frequency resolution which is not always desired.
A wavelet transform is an alternative method for identifying the frequency content of a signal where the window size typically varies based on the frequency, thus creating an inverse relation ship between time localization and frequency localization. That is: high-frequencies have good time localization but poor frequency localization, and low-frequencies have good frequency localization but bad time localization.
Let's say I want to use a wavelet transform as an alternative to the STFT. What specific wavelet should I use? Do some wavelets work better for encoding frequency content?


texts: Octave: Poor performance of convolution by separable 1D kernels

texts:
As a newbie, I wrote a custom function in octave to perform a 2D image convolution using separable kernels.
The results of this custom function were compared with conv2() and they were consistent. But that was where the joy ended.

conv2() works at the speed of light while my custom function is slower than a steam locomotive.

How can I speed up the loops below? As you can see, the number of MACs per pixel is 6 compared to 9 that would have been the case in 2D convolution. The function takes as inputs a row vector(3x1) and a column vector(1x3).

function OutputImage = Convolve3X3(InputImage, RowCount, ColCount, Kernel_x, Kernel_y)

% Create a padded image
PaddedImage = uint8(zeros(RowCount + 2, ColCount + 2));

% Create a staging area
StagingImage = uint8(zeros(RowCount + 2, ColCount + 2));

% Create a Row Vector
RowVector = uint8(zeros(1,3));

% Create a Col Vector
ColVector = uint8(zeros(3,1));

% Copy the input image into the padded image
PaddedImage(2:RowCount + 1, 2:ColCount + 1) = InputImage;

% 1D convolution of necessary rows with Kernel_x
for i = 2: RowCount + 1
  for j = 1:ColCount
   RowVector = PaddedImage(i,j:j+2) .* Kernel_x;
   StagingImage(i,j) = RowVector(1) + RowVector(2) + RowVector(3);
  end
end

% 1D convolution of necessary columns with Kernel_y
for i = 1: RowCount
  for j = 1:ColCount
   ColVector = StagingImage(i:i+2,j) .* Kernel_y;
   OutputImage(i,j)  = ColVector(1) + ColVector(2) + ColVector(3);
  end
end

endfunction

How is conv2() so radically fast?


texts: Does a 20 kHz sine wave sampled at 44.1 kHz need to be reconstituted before playback?

texts:
I am interested in building a synthesizer but I am confused about the Nyquist-Shannon Theorem and sampling rates.

If I sample a 20 kHz sine wave at 200 kHz and plot the data points, the plot looks similar to that of the true signal. But, many websites I have read reference the Nyquist-Shannon theorem and say audio only needs to be recorded at twice the frequency or higher of the highest frequency one expects to record (or 40+ kHz since humans cannot hear audio frequencies greater than 20 kHz). When I plot the data from sampling the wave at 44.1 kHz, it looks nothing like a sine wave.

Before playing a 20 KHz sine wave sampled at 44.1 kHz on my speakers, must I do something to to reconstitute the original signal from the samples?

What confuses me is that I have seen microcontrollers play audio (.wav files) without reconstructing a signal, simply by taking sample data and playing it at the sample rate -- but I don't think that method will work if the frequency of an audio signal is close to half of the sampling frequency. I have found a lot of literature about sampling rates but nothing about playback rates.


texts: How do I calculate the step response of a discrete-time system?

texts:
I am given the impulse response

$$h[n]=2\left (\frac{1}{2}  \right )^n\sigma [n]-\delta [n]$$

where $\sigma$ is the unit step function and $\delta$ is the unit impulse function. I know that

$$a[n]=\sum_{k=-\infty}^{\infty}h[k]\sigma [n-k]=\sum_{k=-\infty}^{n}h[k]$$

Thus, my step response should be

$$a[n]=\sum_{k=-\infty}^{n} \left( 2 \left (\frac{1}{2}  \right )^k \sigma [k]-\delta [k] \right)$$

How do I calculate this?


texts: Find $v_3$ such that $\{v_0, v_1, v_2, v_3\}$ is an orthogonal basis of $\mathbb{R}^4$?

texts:
In this question, we consider the Hilbert space of vectors in $\mathbb{R}^4$ with one of its basis $\{v_0, v_1, v_2, v_3\}$.
Given the first three basis vectors
$$v_0=\left[\begin{array}{c}
1/2\\
1/2\\
1/2\\
1/2
\end{array}\right]v_1=\left[\begin{array}{c}
1/2\\
1/2\\
-1/2\\
-1/2
\end{array}\right]
v_2=\left[\begin{array}{c}
1/2\\
-1/2\\
1/2\\
-1/2
\end{array}\right]$$

How many possibilit(y/ies) are there for $v_3$, such that $\{v_0, v_1, v_2, v_3\}$ is an orthogonal basis of $\mathbb{R}^4$? Give one example of $v_3$ please.


texts: Optimal $ n $ -th Order IIR /AR Approximation of a Moving Average Filter

texts: I would like to approximate a moving average filter with an IIR filter of much lower order than the tap-length of the moving average filter. Optimality shall refer to the $L_2$ norm of the impulse response (up to changes in the total gain). Ideally, I'd also like to constrain the impulse response to be monotonic (i.e. no ripple), but I fear that this requirement along will already fully constrain the filter. If that is the case, I might be willing to take some ripple in turn for a sharper roll-off.


texts: What are the Elements of Filters for Highly Volatility Price Series with Jumps?

texts:
What are the best filtering techniques for a highly volatility price series with jumps. What needs to be considered in designing new filters for practical finance which are adaptive and highly reactive to jumps?

Also I have seen a filter called JMA which is proprietary. If I want to create a similar filter, what direction should I take?


texts: What is the intuition of "averaging is a low pass filter"?

texts:
As I know, the shape of a low pass filter in time and frequency are as follow:

But how averaging work like a normal function in time domain?
Is there any intuition in it?


texts: What are some ways to detect the fundamental frequency in a signal using without using FFT?

texts: I am trying to estimate the fundamental frequency in a signal to corroborate the result which I have got by FPGA implementation of an FFT algorithm. Hence, could someone give me some idea how to do so?


texts: Why multiplying frequency kernel with sinusoid results shifting frequency domain

texts:
I am having great time while reading "The Scientist and Engineer's Guide to Digital Signal Processing" book.

In "Chapter-14 Introduction to Filters" author indicates we can create high pass filter from low pass filters by multiplying low pass filter's frequency kernel with a sinusodial which has 0.5 frequency. Link to the this chapter

I thought I read every chapter carefully but I can't understand why this will work?
Can somebody please explain the reasoning behind of this process as easy that a computer science graduate could understand.


texts:
I know of some oscilloscopes (DSA8300) that repeatedly sample at a few hundred kS/s to reconstruct a few GHz signal. I was wondering if this could be extended to 2D signals (photographs). Can I take a series (say 4) of still pictures using a commercial 16MP camera to finally reconstruct a 32MP image? Will doing this remove aliasing I have from each image?

Should such a thing be attempted from a single image it would obviously not work as no new information is being introduced. If all the pictures taken are absolutely identical, will I still be at the same point as having one image? So are variations essential? Is CCD / CMOS noise enough variation for such a thing to work?

Is there a name for such a technique or algorithm? What should I look for?


texts: Does this system of equations have a solution?

texts:
My heuristic approaches to find a solution for the following system of equations have failed so far. Does a solution exist?

$$
  \left[ {\begin{array}{cc}
   a_k & 0 \\
   a_{k-1} & a_k \\
   0 & a_{k-1}
  \end{array} } \right]
\left[ {\begin{array}{cc}
   g_{k-1} \\
   g_k
  \end{array} } \right] = \left[ {\begin{array}{cc}
   0.5 \\
   0\\
   0.5
  \end{array} } \right]+\left[ {\begin{array}{cc}
   d_{k-2} \\
   0\\
   d_k
  \end{array} } \right]
$$
such that $\{a_k\}$ are known data symbols with equal probability, chosen from L equidistant values. Unknown $\{g_k\}$ and $\{d_k\}$ are to be found such that $E\{d_k\} = E\{d_{k-2}\}=0$.


texts: Adaptive Band pass filter for extracting respiratory waveform from accelerometer data

texts:
I am recording data from an accelerometer attached to the chest (1000Hz). I need to extract the respiratory waveform. I tried an adaptive bandpass filter based on a dominant frequency in my signal based on https://cutt.ly/j24zdj.

Steps in brief: Take frequency spectra, find dominant freq f0, Make a bandpass filter as [max(0.1,f0-0,4), 0.4+rm]

The respiration rate could be between anything from 0.1 to 2 Hz. For deep and normal breathing, the filter works just fine. But for rapid and shallow breathings (low amplitude but high freq), it fails. It identifies the dominant freq incorrectly and thus the entire filter goes wrong.

For example in the image below: The dominant freq should have been around 0.6, but is identifies as 1.7 (Slight diff in amplitude). So filter works badly.

1) How do I design a filter for this? Any suggestions. I can't control the noise in this freq range. The person is as still as possible.

2) The red line is spectra after filtering. Why has the amplitude of certain frequencies increased? I am using the 'spec' function in R: uses the hanning window for computing.

It's not the issue with plotting. The amplitudes have actually increased.


texts:
so, first let's do this for the circular correlation case (it might be easier).

let's say we have three real and periodic functions or signals $x(t)$, $y(t)$, $z(t)$ all with period $2 \pi$:

$$ x(t+2\pi) = x(t) \quad \forall t $$
$$ y(t+2\pi) = y(t) \quad \forall t $$
$$ z(t+2\pi) = z(t) \quad \forall t $$

define the cross-correlation between any two, say $x(t)$ and $y(t)$ as

$$ R_{xy}(\tau) \triangleq \frac{1}{2\pi} \int\limits_{-\pi}^{\pi} x(t) y(t+\tau) \ dt $$

let's say that the maximum value of $R_{xy}(\tau)$ occurs at $\tau=\tau_{xy}$

$$ R_{xy}(\tau_{xy}) > R_{xy}(\tau) \quad \forall \tau \in (\tau_{xy},\tau_{xy}+2\pi) $$

(since $x(t)$, $y(t)$, and $z(t)$ are all periodic with period $2 \pi$, so also are any cross-correlation between any pair.)

$$ R_{xy}(\tau_{xy}) = R_{xy}(\tau_{xy} + k \ 2\pi) \quad \forall k \in \mathbb{Z} $$

so here's the question: except for a difference of some multiple of $2\pi$, is

$$ \tau_{xz} = \tau_{xy} + \tau_{yz}$$

or

$$ \tau_{xz} = \tau_{xy} + \tau_{yz} + k \ 2\pi $$ for some integer $k$?

i think this should be true, but it would be nice if it is proved to be so.


texts: Doubt regarding definition of pixel bond

texts:
The following is a quotation from an Image processing textbook:

The connectivity relationship between a center pixel and its eight neighbors can be quantified by the concept of a pixel bond, the sum of the bond weights between the center pixel and each of its neighbors. Each four-connected neighbor has a bond of two, and each eight-connected neighbor has a bond of one.

I am searching different sources to understand what is pixel bond. How four-connected neighbor has two and eight connected has one.

What exactly is the bond weight?

(The passage from this book seems to also be available via Google Books)


texts: Kernel vs impulse response: what's the difference?

texts:
I am reading The Scientist and Engineer's Guide to
Digital Signal Processing by S.W. Smith. In Chapter 9: Applications of the DFT in the 3rd paragraph he writes "kernel (impulse response)"
What does it mean? What is meant by kernel? Is it same as impulse response?


texts:
I've been trying to get DFT from an impulse response, such as $(n+1)[u[n]-u[n-5]]$, but I don't have any clue of where I should start from
Thanks


texts:
I am studying Perfect reconstruction of 2 channel filter bank for using Biorthogonal Spline wavelet for decomposing ECG signal. I am using Wavelets and Filter Banks By Gilbert Strang, Truong Nguyen to study perfect reconstruction. I am trying to understand the remove distortion condition, and it is given as:

\begin{equation}\label{eqn:updated_no_dist}
{P_0(z) - P_0(-z) =  2z^{-k}}
\end{equation}

Now,it states, they need to design a low pass filter which satisfy the above equation, and to design a low-pass filter, they use this equation.

\begin{equation}\label{eqn:updated_no_distortion}
 {P_0(z) = (1 + z^{-1})^{2p}Q(z)}
 \end{equation}
where p can be any integer and $Q(z)$ a polynomial of degree $2p-2$.

They use $p = 2$ and $Q(z) = -1 + 4z^{-1} - z^{-2}$. I am unable to find the reason why they use $p=2$ and where from the equation $Q(z)$ comes and what is the reason to use that exact equation. I read other papers as well, but I didn't find any reason for choosing $p=2$.

The other question is, what is the relationship between perfect reconstruction and Biorthogonal Spline wavelet.

It would be great if someone can explain this.


texts: How to get the conjugate symmetry of the complex matrix using MATLAB?

texts: How to get the conjugate symmetric of the complex matrix using matlab


texts: Confusion in understanding the Proof of DTFT

texts:
While understanding the proof of DTFT from Signals and Systems by Oppenheim, I have confusion in understanding few steps.

$$ x'[n]=\sum_ {k=<N>} a_ke^{jk(2\pi/N)n}$$
$$ a_k= \frac{1}{N} \sum_ {n=<N>} x'[n]e^{-jk(2\pi/N)n}$$
Since $x[n]=x'[n]$ over a perid that includes the interval $-N_1 <= n <= N_2$, it is convinent to choose the interval of summation to include this inteval, so that $x'[n]$ can replaced by $x[n]$ in the summation. Therefore,
$$ a_k= \frac{1}{N} \sum_ {n=N_1}^{N_2} x[n]e^{-jk(2\pi/N)n} = \frac{1}{N} \sum_ {n=-\infty} ^{+\infty} x[n]e^{-jk(2\pi/N)n}$$
While in the second equality we have used the fact that $x[n]$ is zero outside the interval $-N_1 <= n <= N_2$

I have confusion in understanding that how summation from $-N_1  N_2$ changes to $-\infty  +\infty$. The statement they mentioned above is so confusioning for me.


texts:
The convolution Theorem of the Hartley Transform is given by

$$
Y_{Ha}(k) = X_{Ha}(k)H_{Ha}^{e}(k) + X_{Ha}(-k)H_{Ha}^{o}(k)
$$

where,

$$
H_{Ha}^{e}(k) = \frac{H_{Ha}(k) + H_{Ha}(-k)}{2}
$$
$$
H_{Ha}^{o}(k) = \frac{H_{Ha}(k) - H_{Ha}(-k)}{2}
$$

I tried the following implementation in Python,

HIMG_ = np.flip(np.flip(HIMG, -1), -2)
HKERNEL_ = np.flip(np.flip(HKERNEL, -1), -2)

HKERNEL_e = 0.5 * (HKERNEL + HKERNEL_)
HKERNEL_o = 0.5 * (HKERNEL - HKERNEL_)

ya_hartley_1 = HT((HIMG * HKERNEL_e) + (HIMG_ * HKERNEL_o))
imshow(ya_hartley_1)

And this was the output,

This is the output if convolution is done on spatial or Fourier domain,

The complete code is available on a kernel on Kaggle

Here is my question,

  What is wrong with it? Why the output is not the same as the spatial or FFT convolution?


texts:
I have a time series with lots of steps/jumps (data file here). A plot is given below. I would like to subtract an appropriate value for each of these square wave features to bring them back down to the baseline of the signal. A median filter works really well for removing a small number of outliers in a row, but in this case I probably need a different approach since the square wave jumps can have different durations as seen. A common method I've seen for doing this is to compute first differences of adjacent samples, and look for large differences to detect jumps. I implemented this method but the problem is it often fails, since the one tunable parameter for the method is a threshold value $t$ which the first differences must cross in order to detect a jump:
$$
| x_{i+1} - x_i | > t
$$
As can be seen in the plot below, the jumps I have are often different sizes, so a constant threshold value isn't the best approach. In particular, in some cases there is an interesting signal where adjacent samples can change by large values without being a jump! I have highlighted such a region in red.

Below is a zoomed in view of the red box area. You can see there is a square wave jump followed by an interesting signal. The red arrows depict a place where adjacent samples from an interesting signal have a larger distance between them than some of the jumps in the signal. Therefore a constant threshold method with finite differencing will not work for me.

Does anyone know of a robust procedure to detect and subtract the square wave jumps to end up with a smoothly varying signal with no jumps? I'm sure this must be a solved problem but I haven't had much luck searching online.


texts:
I'm trying to simulate an OFDM signal that has been captured with a receiver that has a much larger bandwidth than the OFDM signal.  I tried the following MATLAB code to simulate an OFDM signal that's 20% of the captured bandwidth and then demodulate it:

Nsymbol = 100;
M = 16;
Nfft = 128;
CPlen = 12;
% random bits...
in_bits = randi([0 1], log2(M), Nsymbol*Nfft);
% use a 16-QAM symbol constellation
symbols = qammod(in_bits, M, 'InputType', 'bit');
% reshape it ofdmmod
parallel_symbols = reshape(symbols, Nfft, Nsymbol);
baseband = ofdmmod(parallel_symbols,Nfft,CPlen);
% oversample by a factor of 5 and then decimate by a factor of 5 again
oversampled = resample(baseband, 5, 1);
back2baseband = resample(oversampled, 1, 5);
symbols2 = ofdmdemod(back2baseband, Nfft, CPlen);
plot(symbols2(:), '.');

When I run the code, I expect that the symbol constellation should be (roughly) the same as the original 16-QAM constellation, but instead I get this:

As a result, I get a BER of about 1% with no noise in my simulation... clearly I am doing something wrong.  What is the right way to demodulate an oversampled OFDM signal in MATLAB?  I tried using the DigitalUpConverter and DigitalDownConverter objects as well, but that was a total disaster.

EDIT: Thanks to Dan's answer below, I realized that the filters that MATLAB applies when downsampling induced a linear phase shift so I was looking at the wrong sample.  When I just do back2baseband = oversampled(1:5:end) I recover the ideal constellation.


texts: Gain at given frequency from $z$-plane zero-pole plot. Two methods gives different results

texts:
I have two zeros at $z=-1$ and two complex conjugate poles at $z=A\cos\theta\pm jA\sin\theta$

This gives me the next transfer function

$$H(z)=\frac{1+2z^{-1}+z^{-2}}{1-2A\cos\theta z^{-1}+A^2z^{-2}}$$

To get filter's DC gain I'm trying to substitute $z$ with $1$:

$$G=\frac{4}{1-2A\cos\theta+A^2}$$

Following graphical method described here (gain as a product of lengths of vectors from every zero to given point divided by product of lengths of vectors from every pole to given point) I've got next formula:

$$G=\frac{2}{\sqrt{1-2A\cos\theta+A^2}}$$

Where is my mistake? Why I've lost square root in the first case?


texts: FIR Filter Design - input parameters conversion

texts:
I am new to DSP and filter design. I have developed a code in C++ to calculate FIR coefficients using Parks-McClellan algorithm. The inputs to calculations are:

Filter type (Low-Pass, High-Pass)
Passband Freq [Hz]
Stopband Freq [Hz]
Passband Ripple [dB]
Stopband Attenuation [dB].

I got a new requirement for inputs, namely to use:
1. Cutoff freq (6dB@1.6kHz)
2. Slope [24dB/octave].

Is there a way to convert these inputs into the inputs I used in my C++ code, as listed above?

Or, more general question, how to design FIR filter with custom slope?

Thank you.


texts:
I need to calculate frequency response, phase response and apply to signals the Butterworth, Chebychev1 and Chebychev2 band-pass filters.

I'm developing in C++ with Qt, and I'm looking for algorithms or a (free) library to do that. I have my plotting solution, I just need the formulas.

The input values I have are the time/amplitude signal, the order of the filter and the low and high cut-off frequencies.


texts: Phase shift of discrete cosine transform (DCT)

texts:
The most common type of discrete cosine transform (DCT-II) is defined as
\begin{align}
X_k&=\sum_{n=0}^{N-1}x_n\cdot \cos\left(\frac{\pi}{N}\left(n+\frac{1}{2}\right)\cdot k\right)&\text{where }& k=0,1,...,N-1
\end{align}
(c. f. https://en.wikipedia.org/wiki/Discrete_cosine_transform)
My problem is that I just don't get what that phase shift of $\frac{1}{2}$ is for. In the paper introducing the DCT-II it was just defined that way, seemingly because the authors wanted to link it to Chebyshev polynomials. That analogy may be neat, but what are the actual benefits of the phase shift?
For all I can say, we might as well leave it out. But obviously, there must be something to it since it is still used today. Sadly all the books I've consulted just give the definition of the DCT without justifying the phase shift. Could someone explain or point me to a good source?


texts: Signal decomposition: getting phase with known amplitude and frequency of main components

texts:
Suppose you have sampled a signal and you know amplitude $ \alpha_i $ and frequency $ f_i $ of the main components for the decomposition:

$$
y(t) = \sum_{i=0}^N{\alpha_i\cos(2\pi f_i t + \phi_i)}
$$

Which methods would be the fastest and the most accurate to calculate the phase vector $ \phi_i $?

Boundaries:

Signal sample has constant sampling rate
Signal sample length is not guaranteed to be an integer multiple of the signal period, but is longer than one period

EDIT

The simples way I can think of is to minimise the RMS error with some optimisation algorithm as shown in code below. I’m wondering if there is a more efficient and simpler way to achieve the same.

%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as op

# Build a signal as sum of cosines
def f_signal(t, magnitude=[], frequency=[], phase=[]):
    y = 0
    for m, f, p in zip(magnitude, frequency, phase):
        y += m*np.cos(2*np.pi*f*t+p/180.*np.pi)
    return y

# Length [s] and number of samples of the signal
t_length, N = 0.09, 1000

# Equally spaced time axis
t = np.linspace(0, t_length, N)

# This is the reference signal (we know all information)
y = f_signal(t, [8, 4], [80, 240], [-90, 40])

# Now suppose we lost the phase information and want to reconstruct it

# Lets use the root mean squared error as function to be minimised
def f_rms_error(phase):
    return np.sqrt(np.sum((y - f_signal(t, [8, 4], [80, 240], phase))**2)/len(t))

# Initial guess
phase0 = [0, 0]

# Minimize RMS error with fmin
phase1 = op.fmin(f_rms_error, phase0)

# Print result
print(phase1)

# Reconstruct signal
y_r = f(t, phase=phase1)

# Plot results
plt.figure(1, figsize=[15, 5])
plt.subplot(1,1,1)
plt.plot(t,y)
plt.plot(t,y_r, 'r-‘)


texts:
The amplitude of the signal in Rician fading case is characterized by the Rician distribution (let's call this model 1).  However, the received power in Rician fading is sometimes modelled as (model 2):
\begin{equation}
\label{model2}
P = q + (1-q)p,
\end{equation}
where $p$ is exponential of mean $1$ ,and $q \in [0,1]$ is the part of the energy received in the line-of-sight.

Probability density function for $P$ can be derived as a convolution of $ p \mapsto \delta(p-q)$ and  $p \mapsto u(p)\frac{e^{-\frac{ p}{1-q}}}{1-q}:$
$$\text{PDF}_{\text{model 2}}(p) = \frac{e^{-\frac{ p- q}{1-q}}}{1-q}u(p-q),
$$
where $u(\cdot)$ is a step-function.
 (This can be seen from the fact that Laplace transform of $P$ is a product of Laplace transforms of these functions.)

In the plot that follows, for model 1 we acquire PDF for received power by simulating the Rician distribution for amplitude and taking the second exponent. We use Rician parameters s.t. $K= 1$, $\Omega = 1$.

The distribution of received power is clearly different for $q \neq 0$ in model 1 and model 2. Case $q=0$ is, of course, the Rayleigh-fading case.

My question is: How do these two models relate? It seems that model 2 neglects the change that signals received from non-line-of-sight paths would add destructively to the line-of-sight signal. Signal power never goes under $q$ (which is the line-of-sight element). In model 1 the signal-power can go to zero.


texts:
Kalman filter achieves convergence of state vector by using sensor observations.

Assuming a sensor such like velocity sensor, giving two axis velocity information in X-axis as well as Y-axis(see edit for one more example to capture essence of my question)

I can always design kalman filter in such a way that X-axis and Y-axis velocity are taken as two different scalar sensor readings and also a kalman filter where I take them as a vector sensor reading.

From both theoretical and practical aspects Is there any thing I am losing by taking X-axis and Y-axis as independent readings?

Or taking them together helps in converging faster?
What is a better approach? And Why?

Edit: Assume one is in a airplane and a camera captures an image and tells latitude and longitude of nadir. With this data I would like to correct my State vector. Now if I use latitude as a different reading and longitude as a different reading than using both of the data. Will I lose anything? As pointed by an answer below if the signal model has an interaction among the components of reading I lose that interaction. Now do latitude and longitude have an interaction? How do I capture this interaction mathematically with H matrix?


texts:
My situation is as follows: i am trying to generate a waveform the hard way, by constructing the samples one by one and then saving the result to a .wav file using Python.

When the frequency is constant, everything is fine: i use $y(t) = \sin(2 \pi \cdot f \cdot t)$.
However, if i change the frequency to be a function of time, things go wrong. If the function is a linear one of the form $f(t) = a + bt$, it still works. But if i choose, for example, $f(t) = 40 + 10 \sin(t)$, the max frequency increases over time, reaching a maximum higher than the expect 50Hz.

I have read something about instantaneous frequency, namely, this: Why does a wave continuously decreasing in frequency start increasing its frequency past the half of its length?. But doing the integral evaluation makes the sound even weirder. And the method i currently have works for linear function of time, so i figured there must be something else wrong.

I also tried to generate chunks of sound in the frequency i need, at each time, and then glue them together. I calculated the period of the oscillation, so that a chunk has as many samples as necessary to make a whole period on that frequency, so that there are no "jumps" between different frequencies. But generates a cracking sound in the sample.

Here is an example of the kind of a JavaScript implementation of the kind of sound i need: http://jsfiddle.net/m7US6/4/.

Any help would be appreciated. Thanks.


texts:
Weighed overlap add Filterbank structure is often used in speech coding.
(https://www.onsemi.com/site/pdf/Filtering_Hearing_Ai.pdf)

I have basic question with respect to analysis.

what is the difference conceptually between:

a) direct analysis: e.g La=R=N1=256  This includes:

extract R=La samples of frame from speech
window it with window of length R=N1
compute its R=N1, point FFT

vs
b) WOLA( let La=256 = 8*32=>k=8,N1=32, R=16    . the diagram for this is:

 1.Extract R samples of speech( R'<'N1)
 2. Put these into FIFO of length La= k*N1
 3. window this by window of length La samples
 4. divide the buffer into k buffers, each of length N1
 5. add the buffers
 6. compute   (block_no)mod(R)  circular shift of buffer
 7. compute N1 point FFT

If the input was a delta function, after analysis, what would be the difference in frequency response of the system(or the window) and what advantage we get by doing WOLA?

thanks  for any help that i can get
sedy


texts: How to measure the "correctness" of a pink noise generator?

texts:
I've implemented a pink noise generator for audio applications. How can I measure how close of an ideal pink noise is my output?

I do not have access to Matlab or Mathematica. I would prefer a solution involving open-source software such as GNU Octave or SciPy.

Related to my previous question: Understanding Voss-McCartney pink noise generation algorithm


texts: Segmenting page with dense input and disconnected separators

texts:
I would like to segment pages of a historic book, which is printed densely and shows disconnected black column separators, see below.

What I would like to arrive at is the following segmentation.

I have been playing with ocropy so far, which makes use of scipy/ndimage. But the tool used is not important, I'm interested in methods (morphology, ...) that can solve this.

Any ideas?