texts: Moving Average Filter's Magnitude Response

texts:
I am working on filtering of data using the Moving average filter in MATLAB.
My data is recorded temperature values on the Y-axis and Time(in seconds) on the X-axis
(I'm getting the magnitude response curve of the Moving average filter as attached)

The shape of the curve seems fair for the moving average filter, according to the theory but,
I want to ask, why am I getting all the magnitude value for MA filter positive? (is this correct?)
Implementing other filters like Butterworth etc. on the same data gives the magnitude values to be negative (image attached)

Here is my MATLAB Code for MA (if it helps)
%T.Temp is table column containing Temperature values in double format.

MovAvg = dsp.MovingAverage(5)

movAvgOut = MovAvg(T.Temp)

fvtool(movAvgOut)

%Gives output the first image I attached

Any help is appreciated.
Thank you.
Ps - DSP is not my field. This work is part of a bigger project i'm involved in, and need some help. Please bear with me if I sound irrational.


texts: Audio transfer function modelling with pink noise + signal

texts:
I'm trying to model a transfer function for a noisy audio system, specifically to measure delayed system response. Before I can confidently apply control, I need to verify that I can exert control over the system by broadcasting a signal and listening for the delayed (and potentially noisy) feedback.

I suspect I can add or subtract a unique signal to pink noise to simultaneously calibrate my transfer function as well as verify that the system is responding to my input.

Is there any prior art on this? Is pink noise the right approach here?

Current Plan:

Generate a second of pink noise and calculate its FFT.
Generate a second of unique signal and calculate its FFT. Maybe this can be a simple square wave at 5kHz.
Sum the noise and carrier signal together as test signal. Split the signal into five 200ms chunks, calculate the FFT of the noise
Record a second of (noisy) silence.
Start broadcasting the combined test signal, while still recording.
Take the FFTs of the "empty" and "full" signals. The difference should be the ambient noise.
Look for the source spectrum distribution.
Given that we found the distribution we were looking for, cross-correlate the recording with the source signal to find the delay at k samples.

Does this sound about right? When I'm done, I should know the average ambient noise, the system delay and the harmonic distortion (maybe I need to sweep frequency for this?).

For calibration, performance doesn't have to be real-time. Is there a more complicated sweep function I can use to find a more accurate delay and response?


texts: Is There Any Filtering Technique with Element Wise Product Only Instead of Convolution?

texts:
Generally we apply a filter to reduce noise from signal. For example if the received signal is $S$ we would get a signal $Y$ with reduced noise applying a filter $W$ with the operation Y= convolution(S,W). In other words filter is nothing but convolution. Is there a known method or example of getting a useful $Y$ such that $Y=S.*W$ where $.*$ means element wise multiplication. An example of element wise multiplication is demodulation.

Can anyone give an insight if element wise multiplication can be used to reduce noise under some circumstances?

Edit: Yes, product in time domain means convolution in frequency domain. Still the question under what situation convolution in frequency domain or product in time domain can reduce noise? Is there any such application? Did anyone used such a method?

Edit: There may be a situation where it is much easier to find a vector $W$ such that the element wise multiplication does reduce noise. It may difficult to find the corresponding filter.


texts: Differences in moving window rms due to sampling rate

texts:
I have a 10 second recorded signal band pass filtered from 100-500 Hz. The original sampling rate is 10 kS. Now I run a moving window RMS with a window length: 10 ms on the signal (10/1000 * sampling rate). I then find the samples of the RMS signal that are greater than 2 SD from mean. I get a total number of, say, 4894 samples that satisfy this condition.

Now I down sample the signal to 2 kS, and do the same process with the recalculated window length also of 10 ms (but different sample length due to the sampling rate).
When I apply the same >2 SD criteria I get 974 samples.

I then use the sample indices from the two different cases to see if they are indeed the same samples by normalizing the indices by sampling rate times 1000 to get ms precision.

Then I compute the intersection of these two sets of times when s.d > 2 for both the sampling rates.

I get a value of around 700. I am trying to understand why. Shouldn't all the samples times in the set corresponding to the 2 kS sampling rate actually be a subset of the sampling times seen in the 5 kS sampling rate for the > 2 SD criteria?


texts: Calculate Min/Max value of multiple sine waves added together

texts:
So I feel like this should be possible, I'm just not sure how to do it. Let's say I'm generating a signal with a sample rate of 10K samples/second. Then let's say my signal is comprised of X number of frequencies( we'll use 3 for this example). For each sample I want to take the three frequencies and add them together, lets use 1000Hz, 1250Hz, and 1500Hz. Now that i've added them together, they might be too large (lets say the max value i can output is a 1) so then i need to normalize them. Is there a way to calculate the maximum or minimum value that these three frequencies will produce over a given period of time/number of samples so that i can get normalize my signal accordingly?

Thanks for any input on this!


texts: Amplitude Spectrum, Nyquist Frequency, mixed/min/max wavelets

texts:
The problem is here. Now I know the definition of mixed/max/min phase wavelets, whether the roots lie within the unit circle or not. Starting from n = 1, let
$$ x_t = ( 5, 6) $$
$$ X(z) = 5 + 6z $$
This would have root -5/6, making it a max phase wavelet (since it's within unit circle). If we switch 5 and 6 around in xt, we would get root -6/5, which would make it a min phase wavelet.

I don't know how I would show they have same amplitude spectrum. I know |X(f)| amplitude spectrum would lie between -fc and fc (Nyquist frequency) and fc = 1/2*delta t. Not sure how I would go about proving it.


texts:
I am stuck at question number 2.8

This is how I have gone about solving it:

I have calculated $y(t)$ by convolving $x(t)$ with $h(t)$ using the fact that $x(t)$ convolved with an impulse at $t=t_0$ is just the same signal $x(t)$ time shifted by $t_0$. Hence,

$$y(t)= x(t+2) + 2x(t+1)$$

Now, according to the answer in the book, signal $y(t)$ is $t+3$ on $-2 < t \leq -1$. But, according to my answer, if we add $x(t+2) + 2x(t+1)$, the value of $y(t)$ at $t=-1$ is $4$. Why Please help me.

Also, why is it that in the answer in the book they haven't included $t=-2$ in the range?


texts:
I am currently looking at the STFT for Librosa. I was wondering how to understand the output of the STFT function, specifically what kind of frequencies the different values represent.

Say I have a n_fft of 256, that means that the shape output of the STFT will be 129 (1+n_fft/2). I therefore understand that for each frame I have 129 bins, these bins represent some sort of frequency within this frame. I am assuming that the bins start with a frequency of 0hz and go up to some value, with each bin representing an equally large range.

But how do I figure out what the maximum frequency used for STFT is? And what about the range of each of the 129 bins?

Just to be clear, I have tried looking at the source-code and at the documentation, but I have not really become much wiser.


texts: Effect of Windowing Function on Autocorrelation Function

texts:
I need to compute FFT components and autocorrelation of a discrete signal. When I do not use any window function before taking FFT of the signal I get a reasonable AC function from FFT that includes some information about period of the signal:

However when I apply a window function (Hanning) before applying FFT I got an AC function that seems to be distorted:

What is the effect of window function on an AC function that is computed from FFT?


texts:
I have some basic doubt in understanding a matched filter (MF). Please note that I have already seen very informative post in this context Understanding the Matched Filter

In what follows, I list three information about matched filter that I found in standard textbook and wikipedia:

(1) Wiki page says that MF is obtained by correlating a known signal, or template, with an unknown signal to detect the presence of the template in the unknown signal.
In digital communication, I assume that the known template is the pulse shape p(t) used by the transmitter [e.g., 'sinc-like' impulse response of root-raised cosine (RRC) spectrum] and the unknown signal is the noisy received signal r(t).
Assuming that the receiver knows p(t) used by the transmitter, the receiver convolves r(t) with p*(-t). The output of convolution is sampled at decision points.

(2) A standard textbook on communication says that the receiver does a convolution of the received waveform with the time-reverse conjugate of channel response h(t). That means, MF operation consists of convolution of r(t) with h*(-t), followed by a sampler.
(3) The textbook also says that MF can be seen as an equalizer.

I understand that MF does a convolution followed by a sampling. That is it. Simple. But I am seeking detailed clarity in my understanding of the philosophy of this operation.

Questions:

(A). What is p(t) in the MF operation? Suppose I have NRZ data, then do I use the 'rectangular pulse of NRZ' or the 'sinc-like impulse response of RRC' as p(t) in the convolution? Wiki says, rectangular NRZ pulse is used. Why so?

(B). Refering to the two different notions in (1) and (2), seems my understanding of the relation between p(t) and h(t) is unclear. Can someone make it clear with simple explanation?

(C). How do you explain MF operation as an equalization? Equalizer removes ISI, MF removes the effect of noise...

(D). What is partial MF?

I really appreciate your help in improving my understanding of this beautiful concept of matched filtering.


texts: Advantages and disadvantages of Hilbert-Huang Transform versus Fast Fourier Transform?

texts: Which one is better for frequency detection? In what circumstances? Advantages and disadvantages of both? Which would be better for a Voice Activation Detection algorithm? Which would be better for an instrument tuner? How would one detect pitch user an HHT transform? Is it faster than FFT? How would one implement HHT?


texts: How is time used in this autocorrelation expression?

texts:
There are two papers I'm working through that explain applications of autocorrelation.
The first - by Monti - is quite clear to me:

It runs a sum over the time axis.
The second one by Brown makes no sense to me:

How is time $t$ collapsed as a dimension? Is there an implicit sum? As it stands, the only way I can interpret this is as a multi-dimensional surface with two independent variables, $t$ and $\tau$.


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

texts: How to compensate for the frequency offset of a PAM, PSK, or QAM signal when we are using single carrier transmission, I am looking for algorithms, explanations to be able to carry out the implementation on matlab


texts: Why are polynomial filters for gradients in the opposite sign?

texts:
I was looking at the filters, and the ones for 1st derivatives, gradient http://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter#Tables_of_selected_convolution_coefficients
by testing the 9-points one:

D=[-4;-3;-2;-1;0;1;2;3;4]/60;
y=sin((0:0.2:6)');
yd = conv(D,y);
plot([repmat(0,128,1) y yd(5:end-4)])

we see that actually [4;3;2;1;0;-1;-2;-3;-4]/60 would give the right signs, so why is it like that?


texts: The signal is the white dot or it is a point target. and I have to find the SNR of this image/point target. As i increase noise in thie image, black back ground shows some disturbance. AnHow should i approach this problem. I have to do it in matlab code.


texts: Amplitude of a signal. [-1, 1] vs variants and MFCC implications

texts:
I've seen two sort of audio signals, one where the y-axis takes values between -1 and +1. I'd think this means 1 means loudest? Not really sure

And the second one I saw goes between $-10 000$ and $+10 000$. What do these values actually mean? Are they arbitrary?

I'm doing work on audio features using MFCC but it seems that when I take a file with values between $[-10k,+10k]$, I get normal cepstral coeffecients of which we take the first $12$. But when I take audio with values in $[-1, 1]$, the cepstral coeffecients are in a negative range.


texts: Difference between filtering using Matlab and filtering using C code generated by Matlab

texts:
hi，has anyone used Matlab coder to generate C code for a filter? Here is the thing: I tried to filter the same audio signal using Hamming window-based FIR bandpass filter, the filter is designed by Matlab DSP tool box. Then I used Matlab coder to generate the corresponding C code for deploying the filter on the MCU. But I found that the filtering in Matlab and that in MCU has differences, I think the result from Matlab is correct.
Is the filter C code generated by Matlab not reliable? Maybe I should try ARM DSP library


texts: Relationship between wavelet shape and wavelet filters

texts:
MATLAB has a library of wavelet functions, showing their "continuous forms" as well as the the decomposition and reconstruction filters.

In decimated wavelet transform the filter size remains the same and the data points are downsampled by $2^j$ at every level $j$. A question that I cannot find an answer is how are the wavelet filters and the actual continuous forms are related? I am chemist, exploring wavelets for some applications.

For example, MATLAB has db9. Its continuous estimation, wavelet function, psi or phi do not look similar to the decomposition filter points. One text says that some wavelets start out as filter and later their continuous forms are estimated, but why in the case of db9 and many other, the filter points do not match the shape of the wavelet function?

Thanks.


texts:
We have unnormalized Haar matrix which, is for example,

H4=[1 1 1 1;1 1 -1 -1;1 -1 0 0;0 0 1 -1]

After normalizing it, we use for haar transform. I know how to produce it using kronecker product. Can we generate other unnormalized wavelet family matrixes? For example, how can I generate unnomalized Daubechies or symlet matrixes?

Thanks


texts:
I'm just learning about wavelets and the PyWavelets package.  I saw a reference to MODWT, which led me to the SO post here and then to this Python package: https://github.com/pistonly/modwtpy
What seems to be unique about the modwt and modwtmra functions (compared to anything I'm seeing in PyWavelets) is that their outputs can be summed across the zero axis to reconstruct the original series.
Here's an example:
import numpy as np
import modwt

x = np.arange(8).astype(float) # array([0., 1., 2., 3., 4., 5., 6., 7.])
wavelet = "db2"
wt = modwt.modwt(x, wavelet, level=3)
wtmra = modwt.modwtmra(wt, wavelet)
x_recon = wtmra.sum(axis=0).round(5) # arbitrary rounding to ensure precision of assert
assert (x == x_recon).all()

If I understand correctly, an advantage of this is that the outputs of modwtmra can be used for regression purposes.  That is, the individual output series can be used as model inputs and then the outputs of each model (one for each series) can be summed to make a prediction on the original series.
I see in the PyWavelets documentation that the Stationary Wavelet Transform is similar to MODWT if norm=True, however, I'm not sure if it's possible to get the same kind of additive output.
So what I'm wondering is:

If my understanding about using the outputs from MODWT/MODWTMRA for regression is correct.  If so, could someone explain how exactly this works?  What are these transforms doing that allows the wavelet properties to be maintained while still allowing for regression?
Is it possible to replicate the additive output using SWT in PyWavelets?  If not, what is the difference between what SWT is doing in PyWavelets and what's happening with MODWT/MODWTMRA?


texts:
My starting point in what follows is a radially symmetric random field. Taking the Fourier transform of this (and plotting it in logarithm to highlight the patterns), I obtain the following image in Fourier space:

As you can see, there is a radially symmetric part of concentric circles, superimposed with a cross-pattern. Now, I do not understand this last part, but I highly suspect this to be an artifact that is not supposed to be there...

It would not surprise me if this is a problem that more people have run into this problem, but I have not been able to find an answer yet.

So, bottom line: Why is there a cross-pattern in my image?


texts:
I have a square wave, with 10% duty cycle. It is 90ms LOW, 10ms HIGH. Sampling frequency is 1kHz.

In spectral domain, after following this example, I get fundamental frequency (10Hz), as well as its harmonics. Also, the whole spectrum is enveloped with harmonics if 100Hz (which comes from pulse duration 10ms)

That is simply the case of spectral decomposition working properly for square wave, because of Fourier transform definition, right? I also read similar question but still a bit confused


texts: Extrinsic Camera Parameters without introducing control points?

texts:
I am new to the field of Camera Calibration. We have successfully calibrated our camera through classic method (using checkerboard etc).
I recently realized there is a field of research for Automatic Camera Calibration. I expect these methods to only be capable of estimation intrinsic camera parameters. However, from my early researches, it seems to me, that they can estimation extrinsic camera parameters (X,Y,Z position in the world coordinates). It occurred to me that:

If these methods does not need any control points (real world points with measured distances)?
If it is so, how they do that, I mean how they know where is the camera in the world?


texts: Is there a general name for this kind of filter?

texts:
I have a signal with some high frequency noise.

At first was just adding a low pass filter which removes the noise very well but has two drawbacks: First there is a delay (which I try to avoid in my application) and second, if the input signal changes suddenly (from 10% to e.g. 50% of the maximum input level), the filter removes the high frequency part too and the transition is smoothed out and/or overshoots. This is both expected.

A second idea I implemented is just a kind of bit reduction (in my case a fraction of a few bits). It seems that currently this is exactly what I need. Noise is perfectly removed, there is zero delay, and high frequencies are still preserved.

My only question is, is there a name for this kind of filter (to investigate some further ideas)?


texts:
I have to find The energy of a signal using Rayleigh th. the signal is $$ x(t) = A e^{-At } u(t) $$ assuming A>0

Using the classic definition of E , I found that it should be $$ \frac{A}{2} $$

Using Rayleigh I should do
$$  \int_{-inf}^{inf} | \frac{A}{A+i2\pi f}|^2 df $$
This because I previously found X(f)
From this I obtained
$$\frac{A^2}{8 \pi^2  f}  log (\sqrt (A^2 +4\pi^2f^2)) $$ with log from -infinite to infinite.

$$ log \sqrt x = log \frac{x}{2} $$ but now I have no idea how to continue , and if this is correct.. thank you


texts: Relationship between windowed fourier transform and wavelet transform

texts: I was reading on windowed fourier transform and wavelet transform, and i was thinking that the windowed fourier transform is a subset of wavelet transform. Is that true?


texts:
I have an asynchronous multi-user communication with $10$ users. With a matched filter, I see the average BER by user is small ($<0.08$).

Do you know any methods that can reduce the BER by user on reception? I know the multi-user detection but I don't want to estimate the different channels.


texts:
In this paper decorrelating detector we have the decorrelating detector. With $K$ users, the correlation matrix of the signature waveform is $R=S^TS$ with $S=[s_1,...,s_K]$. Moreover $size(R)=K \times K$

This is my question: In the article (p $679$), we have: $d_1=\sum_{i=1}^K [R^{-1}]_{1,i}s_i$. But, how can we make this product, the elements do not have the same dimensions? In fact, $size(R^{-1})=K \times K$ and $size(s_i)=N \times 1$ with $N \neq K$


texts: FMCW Radar Signal Processing Flowchart / Ambiguity Function

texts:
I Read many papers on signal processing for pulsed and FMCW radars, taking the FMCW type, the signal processing start after mixing the received and transmitted signal, the result of this operation is a signal called FMCW beat signal. We apply 2D-FFT to extract targets from the beat signal then we apply a detection algorithm to recognize true targets and reject others.
I read also from some papers that they use a tool called ambiguity function, I don't know where this function is situated in fmcw signal processing flowchart or this function is just with pulsed radar!
My question is :
What's the role of ambiguity function in radar signal processing? and when we should use it?


texts: Why Z-transform is considered as separate transform?

texts:
The mathematical formula of the Laplace and Z transforms are same with just one difference.  I.e. in the first we use $t$ for continuous-time signal and in the latter uses $n$ for discrete-time signal. I don't think that there are any other differences.

While discussing Fourier transform for continuous-time (normally analog) signals, we use the continuous Fourier transform and for discrete-time (most often digital) signals we use the discrete-time Fourier transform (DTFT).

So my question is why the Z transform is considered a separate transform? Why it is not named the "Discrete Laplace Transform" ?


texts: Fourier transform of the magnitude of the fourier transformed signal

texts:
I've come to realize that the Fourier transform of an already Fourier transformed signal gives the time-reversal signal.
$$\mathcal F(\mathcal F(x(t)))=x(−t)$$
ref 1,ref 2
However, my question is, if Fourier transformation is applied to the magnitude of the Fourier transformed signal (excluding phase information), is there a way to interpret the returned signal?
$$\mathcal F\{|\mathcal F\{x\}|\}=?$$


texts: Time domain maximum from frequency domain data?

texts: Is it possible to calculate the maximum value of a time-domain signal from frequency-domain representation without performing an inverse transform?


texts:
I have a couple of questions regarding windowed FFTs:

Why is the noise floor higher with windowed FFTs (according to Wikipedia's spectral leakage page, anyway), when the whole point of windowing is to reduce side lobes?
I realize that different windows are better for different things, but is there a window that is considered to be the best all-around window for spectrum analysis?  Alternatively, is there a better way to do spectrum analysis than windowed FFTs?  It would have to be a DSP approach (i.e. I can't do an array of analog filters), but within that constraint I am game for different solutions.

Thanks for your time in reading this.


texts:
I am trying to understand the summary of Sift Keys presented here.
From what I understand, we are concerned with finding extrema in scale-space of a difference of gaussians pyramid. So we have an image with some resolution say $(H,W)$, and we repeatedly smooth and sub-sample the image separately with two gaussian filters where the standard deviation of the filter is different (hence the difference), and then take the difference image. Say we have $k$ total levels.
The algorithm then looks for local extrema both in resolution and in scale space, over a (3,3,3) neighborhood, but I do not understand how you can compare neighborhoods/locations at different resolutions.
For instance if we have say (600, 600), (300, 300), (150, 150) as our image resolutions and there is a local max at (200, 200) in a 3 by 3 neighbourhood in the (600, 600) resolution image, do we then simply map this point to (200/2, 200/2) = (100, 100) in the second image, and (150/4, 150/4) = (37.5, 37.5) in the third image (and from there either take floor(37.5) or do some kind of bi-linear interpolation)?
Any insights appreciated.


texts: FFT radix-2 in audio - resolution issues - logarithmic scale - sample rate VS buffer size

texts:
I would like to ask how to deal with audio DSP. Is FFT good at all to deal with audio? I am working on an app, that generates a single impulse, rout it to an EQ plugin, and the output of the plugin is run through the FFT to see what exact changes were made by the plugin.
And I have found some issue with FFT, that I'm not sure how to solve:
FFT gives you the same resolution for the whole frequency range. If I have a sampling rate of 48 kHz, and a buffer (and DTFT) size of 1024 samples, then the resolution is 48000/1024 = 46.875 Hz. And the same resolution is for lower, medium and highest freq, in one word for the whole range. But in audio, we often work on logarithmic scales. So we need better resolution on low freq, than high freq. How do you deal with that?

Of course I can store data from buffer size until I have 48000 samples, and then run it through the FFT, then I would have resolution 1 Hz. Great, but at first: then I can refresh at least only one time per second, at second: I still think I don't need such a high resolution for high frequencies, and at third it didn't solve my problem if we have other sample rate like 44100 which is not 2^L.

So, is FFT not a good choice for that purpose, please tell me in which situations it is good to use radix-2 FFT in audio industry? As I know for filtering there are also better choices like FIR filters. So what is about all that FFT? :)


texts: Covariance between real and imaginary parts of Fourier transform of a stationary time series

texts:
Since Fourier transform of a random stationary process in time (in the case of existence) is not necessarily real, my question is what is the relation between the covariance of real and imaginary parts of the transformed random variable.

More specifically, according to Wiener-Khinchin theorem for the time function $x(t)$ with some conditions the Fourier transform $X(\omega)$ exists.

The question is can one recover  $$\mathbb{E}\left[\Re\{X(\omega)\},\Re\{X(\omega')\}\right]-\mathbb{E}\left[\Re\{X(\omega)\}\right]\mathbb{E}\left[\Re\{X(\omega')\}\right]$$
and $$\mathbb{E}\left[\Im\{X(\omega)\},\Im\{X(\omega')\}\right]-\mathbb{E}\left[\Im\{X(\omega)\}\right]\mathbb{E}\left[\Im\{X(\omega')\}\right]$$

from $$\mathbb{E}\left[\bar{x}(t),x(s)\right]-\mathbb{E}\left[\bar{x}(t)\right]\mathbb{E}\left[x(s)\right]\;?$$

I would be happy to be guided to a reference as well.

This is the most general case but the simpler case where the input signal is real is the main concern for me. I would be happy to be guided to a reference as well.


texts:
I'm trying to develop a prototype at the moment that includes a sensor to detect how fast someone is typing on a keyboard. Essentially this sensor is a microphone that will be able to detect keyboard taps and then count the number of taps in a certain time interval. So far I've just been using the instantaneous energy of the signal, and then counting the high energy impulses. The issue with this is that every sound is recorded, meaning that someone talking over the keyboard taps results in a false reading.

I need to try and design a filter that will only pass the keyboard taps and block any of the other sounds it might hear. I've read about matched filters, but am struggling on how to actually create one. I know I'll need to convolve some characteristic 'keyboard tap' signal with the microphone recording, which should give me a signal with impulses at the start of each of the taps.

I've been using matlab to try and analyse these signals, but now need to turn that into a matched filter, which has left me pretty stumped! Any help on how to design one and any pointers on convolving signals on Arduino would be really appreciated!!!

Thanks,

Luke

***I'm trying to design an intensity tracker for PC gamers btw


texts: Some doubts in the matlab code of computing Wigner distribution

texts:
I found that a discrete form of the Wigner distribution was
$$WD(n,k) = \sum_{m = -N/2}^{N/2}f(n+m)f^*(n-m)e^{-j\frac{2\pi}{N+1}2mk}$$
where $f(n)$ is the signal and time limited within $|n|\le N/2$. I manage to find the matlab code of the Wigner distribution (Sorry that I am not very familiar with the format to put the codes):

 1. function [WD] = EVD_tfrwv(x0)
 2. X=fft(x);
 3. X=[X(1:N/2+1);zeros(N,1);X(N/2+2:N+1)];
 4. x=2*ifft(X);
 5. x=[zeros(N,1);x;zeros(N,1)];
 6. X=zeros(2*N+1);
 7. for k=1:2*N+1   X(:,k)=x(k+(0:2*N));  end
 8. ww=X.*conj(flipud(X));
 9. WW=fft(ww([N+1:2*N+1,1:N],:));
 10. WW=real(WW([N+2:2*N+1,1:N+1],:));

My question is

I understand line 8 tries to get the correlation $f(n+m)f^*(n-m)$ (probably) for different time $n$, but I get confused when line 9 performs the FFT after rearranging the order of these correlation samples. Could anyone tell me why we have to do this?
Line 10 has similar operation.
Why there is a real operator in line 10? I do not see any real operator in the Wigner distribution formula although I understand Wigner distribution is a real distribution.

Thanks for anybody who tries to help!


texts: Difference between Fine and Coarse Frequency Offset in OFDM

texts:
I need to explain for a project the difference between Fine and Coarse Frequency Offset in OFDM signals. I know what the frequency offset is and what it does to the carrier frequencies but I cannot understand where the terms "coarse" and "fine" refer to. Those terms are also used in GNU Radio in the Schmidl and Cox sync block.

Most of Google results take the fine and coarse frequency offset as granted and there is no specific explanation.

Can someone explain this to me?


texts: Why a pair of complex-conjugate zeros provide a nulling filter? (FIR filter case)

texts:
FIR filter having a pair of complex-conjugate zeros that lie on the unit circle, with zeros of the form:
$$
z_1 = e^{j\omega_i}\qquad\text{and}\qquad z_2 = e^{-j\omega_i}
$$
And transfer function:
$$
H(z) = \left(1 - z_1z^{-1}\right)\left(1 - z_2z^{-1}\right) = 1 - 2\cos\left(\hat{\omega}_i\right)z^{-1} + z^{-2}
$$

Has the following pole-zero map, and magnitude & frequency response graph:

But what is the intuitive idea behind that? And how the locations of zeros will affect the result? Thank you!


texts: Horizontal banding (flickering) due to electronic rolling shutters

texts:
It is a well-known artifact that in CMOS cameras with electronic rolling shutter, horizontal banding (flickering), i.e. brightness intensity variations, are observed when the image is recorded under fluorescent lighting. This is due to sequential exposures of the rows of the image, which leads to different intensities in rows of the image due to the 50 Hz (Europe) or 60 Hz (US) AC power line frequency. The light flicker frequency is twice the AC power line frequency, i.e. 100 Hz or 120 Hz, respectively. The figure below (Yoo, 2014) shows the flicker generation process in a CMOS image sensor:

A common solution proposed in the literature for this artifact is to synchronize the exposure time (shutter speed) to an integer multiple of the light modulation period. For example, as seen in the figure, we can set the shutter speed to 1/60 s, and hence avoid flickering (1/120 s would work also since the light modulation frequency is twice the AC power frequency!). However, this solution requires detecting the presence of flicker and estimating the flicker frequency. Various DFT based methods are proposed in the literature but they are usually based on the difference images produced from sequential frames, and are known to not work very well with single images (still frames captured by a digital camera).

I would like to know if there are any image processing methods which can be applied to a single image, which shows horizontal banding, to remove the horizontal bands through post-processing. In other words, horizontal and periodic intensity variations in the still image needs to be eliminated. A possible scenario is where the shutter speed can take any value (i.e. no sync of exposure time to light modulation frequency), and a still image with horizontal bands has already been recorded. I would like to take this image, apply some kind of filtering (possibly across columns) and get rid of the horizontal bands. In my literature search, I encountered only one such method which suggests using a variation of classical Gaussian filtering (with different Gaussian variances depending on bin intensities), but I could not yet find any other image processing solution. How else can horizontal bands possibly be eliminated through post-processing?

I would appreciate any references to the literature or suggestions for a possible solution.


texts: Problem with the definition of Linearity

texts:
From the high school maths we know that y=mx+c is a linear equation. However, in DSP the linear system must satisfy Additivity properties which y=mx+c does not hold because of +c.
So, is the definition of Linearity different in DSP and Maths? if so why?
Thanks.


texts: How to form Kalman filtering matrices for a problem with variable acceleration?

texts:
Assuming we have

time vector $T$ with constant time step $dt$
position vector $X$
velocity vector $V$
acceleration vector $A$

All vectors $X, V, A$ have noise on their measurement ( $n_x$ , $n_v$ , $n_a$ ). And all should be filtered.

How to form Kalman filter for the the matrices:

$$ x_k= A x_{k-1} + B u_k + w_{k-1} $$
$$ y_k= H x_k + v_k $$

where $w$ and $n$ are normal noise. Each $x_k$ is $\begin{bmatrix} X_k  & V_k & A_k\end{bmatrix}^T$

I came up with two matrix assumption for kalman filtering:

which one is correct?
$$A=\begin{bmatrix} 1 & dt & \frac{1}{2}dt^2 \\ 0 & 1 & dt \\ 0 & 0 & 1 \end{bmatrix}$$
$$B=\begin{bmatrix}0 & 0 & 0\end{bmatrix}^T$$
$$H=I$$

or this one?

$$A=\begin{bmatrix} 1 & dt & \frac{1}{2}dt^2 \\ 0 & 1 & dt \\ 0 & 0 & 0 \end{bmatrix}$$
$$B=\begin{bmatrix} 0 &  0 & 1 \end{bmatrix}^T$$
$$H=I$$

The acceleration is not constant and there is a variable input acceleration. I need to fix acceleration as well as velocity and position.

Please explain why which one is correct and the other is wrong.

This is my code. How to extend it for variable acceleration?

time_resolution=0.016;
tmp1=diag((1:9)*0+1,3);
tmp2=diag((1:9)*0+1,6);
A_T=eye(9)+tmp1(1:9,1:9)*time_resolution+0.5*tmp2(1:9,1:9)*time_resolution*time_resolution;
A=A_T;
C=eye(size(A));
p_tmp=[position_x;position_y;position_z];
v_tmp=[velocity_x;velocity_y;velocity_z];
a_tmp=[acceleration_x;acceleration_y;acceleration_z];
X=[p_tmp;v_tmp;a_tmp];

xh_=X*0;
yh_=X*0;
inov=X*0;
xh=zeros(size(X)+[0,1]);
xh(:,1)=X(:,1);

R=1;
Q=100*diag([.1 .1 .1 .5 .5 .5 10 10 10]);
Px=eye(size(A));

for i = 1 : size(X,2)
    % x(t|t-1) = A*x(t-1|t-1)
    xh_(:,i) = A * xh(:,i);

    % P(t|t-1) = A*P(t-1|t-1)*A' + Q
    Px_ = A*Px*A' + Q;

    % K(t) = P(t|t-1) * C' * inv(C*P(t|t-1) * C' + R) )
    K = Px_ * C' /(C*Px_*C' + R);

    % y(t|t-1) = C*x(t|t-1) + R
    yh_(:,i) = C * xh_(:,i) + R;

    % innovation error  y(t) - y(t|t-1)
    inov(:,i) = X(:,i) - yh_(:,i);

    % x(t|t) = x(t|t-1) + K(t)*(y(t)-y(t|t-1))
    xh(:,i+1) = xh_(:,i) + K * inov(:,i);

    % P(t|t) = (I - K(t)*C) * P(t|t-1)
    Px = Px_ - K*C*Px_;
end


texts: Why only take 205 sample when using a sample rate of 10khz or 8khz?

texts:
I've often seen i my DSP book, that when they record a sound, they use a sample rate at 10khz or 8 khz, which makes sense since the fs is 2 times bigger than the sampling signals max freq.
What i don't understand is why they only take 205 samples, rather than 10000 samples or 8000 samples, wouldn't each sample be "unperiodic"??..
What is the magic behind the number 205 ???

EDIT:
I know they aren't able to sample 8000 times here, but why not take 320 samples then?? than just 205??

Example:

  Estimating DTMF Tones with the Goertzel Algorithm

  The minimum duration of a DTMF signal defined by the ITU standard is
  40 ms. Therefore, there are at most 0.04 x 8000 = 320 samples
  available for estimation and detection. The DTMF decoder needs to
  estimate the frequencies contained in these short signals.

  One common approach to this estimation problem is to compute the
  Discrete-Time Fourier Transform (DFT) samples close to the seven
  fundamental tones. For a DFT-based solution, it has been shown that
  using 205 samples in the frequency domain minimizes the error between
  the original frequencies and the points at which the DFT is estimated.

  Nt = 205; original_f = [lfg(:);hfg(:)]  % Original frequencies
  original_f =

     697
     770
     852
     941
    1209
    1336
    1477

  k = round(original_f/Fs*Nt);  % Indices of the DFT estim_f =
  round(k*Fs/Nt)      % Frequencies at which the DFT is estimated
  estim_f =

     702
     780
     859
     937
    1210
    1327
    1483

  To minimize the error between the original frequencies and the points
  at which the DFT is estimated, we truncate the tones, keeping only 205
  samples or 25.6 ms for further processing.

  tones = tones(1:205,:);


texts: Does there exist a discrete low-pass filter algorithm that is greatly more computationally efficient for high-order (e.g. greater than 10) low-pass filters than straightforward methods like series of IIR Butterworth filters?


texts: Seperate music instruments and estimate note density

texts:
Here is my problem: I have a .wav file of a musical performance with three instruments: a clarinet, a bass and drums. My goal is to detect the note speed played by the clarinet player (so I want the number of notes played every five seconds: ex: between 0s and 5s = 4 notes played; between 5s and 10s = 2 notes played, etc...).
To seperate I used spleeter by Deezer and to calculate the note density, I used the librosa function librosa.onset.onset_detect() however I am note quite satisfied with the result.
Does anybody have a better idea to solve my problem?
Thank you in advance!


texts: What features describe audio signals? (Besides frequency and amplitude)

texts:
I recorded sounds with a microphone and I try to distinguish them in my Java program. The frequency works quite good, but if I look at the fourier transforms it seems like there should be more features to distinguish sounds. I don´t know very much about signal processing, maybe you can help me. Here is a picture of two fourier transformations.

http://i.imgur.com/7npkt.png

I know that the frequency is determined by the index of the maximum magnitude (hope that is the right term). In the first fourier transform curve, its at 100, in the second at 12 (frequencies are 1102.5 and 132.3). But the two sounds look so differently when they are transformed, what else could I use to distinguish them?


texts:
suppose that we have a wave which has some peaks in some indexes like 2,12,25,... and otherwise zero.

I want to extract the time lag or lags between peaks of a signal. could anyone help me ?


texts: Fast method to decompress huffman and arithmetic encoded data

texts:
By nature, when decompressing huffman (and I assume also Arithmetic encoded data) we have to read it in bit by bit and then see how it traverses a huffman code tree that has the actual symbol value at each leaf node. Since processors are designed for byte or multiple byte read/write, this seems to be a very slow and inefficient process.

Since huffman coding and many other methods are used to encode image and video data which produces non-byte boundary aligned results, how do we speed up the decompression process? I would assume that instead of reading 1 bit at a time we read several bits or even bytes and then using some clever method, find all symbols in it at once. I think that this is how all codecs for video and image would work to decompress the compressed image data fast. However, this is only my assumption and I do not know of any details. Could someone kindly reveal the secret to me?


texts: Running Integral of sine and cosine functions

texts:
In typical signal processing course we were taught that the integral of signal $x(t)$ is given by
$$y(t) = \int_{-\infty}^{t}x(\tau) d\tau$$
How can we use this definition to evaluate the integrals of eternal sinusoids, $x(t) = \sin (t)$ or $x(t) =\cos(t)$. I mean how do you evaluate / prove that $$\int_{-\infty}^{t}\cos(\tau)d\tau = \sin(t)$$

Thanks in advance.

Edit: I raised this question because of its (ab)use in frequency modulation(FM) for single tone signal. To give details, any FM wave of single tone signal $m(t) = \cos(\omega_0 t)$ is given by
\begin{align}
\Phi_{FM}(t) & = A\cos\left(\omega_c t + K_f\int_{-\infty}^{t}\cos(\omega_0 \tau)d\tau\right)
\end{align}
It can be written as
\begin{align}
\hat{\Phi}_{FM}(t) & = A\: Re\left\{\exp\left(j\omega_c t + jK_f\int_{-\infty}^{t}\cos(\omega_0 \tau)d\tau\right)\right\}
\end{align}
or
\begin{align}
\hat{\Phi}_{FM}(t) & = A\:Re\left\{\exp\left(j\omega_c t\right) \exp \left(j K_f\int_{-\infty}^{t}\cos(\omega_0 \tau)d\tau\right)\right\}
\end{align}
In communication systems textbooks, the argument of exponential term involving cosine term is magically written as $\dfrac{K_f}{\omega_0} \sin (\omega_0 t)$ which really fascinates me.

Later they claim this is periodic signal and hence can be represented as Fourier series expansion involving Bessel functions etc.,

As I mentioned it is really magic and am interested in knowing this magic trick. Hence this post. Sorry for being so long....

Thanks to all of you.


texts: How to detach overlapping contours but not to remove small ones. Alternative to erosion-dilation cycle

texts:
Background:
I'm trying to detect trees on a sequence of aerial images by calculating a stereo disparity for consequent pairs of images. This is how disparity looks like. White spots are trees:

Next, I am locating contours on the disparity map and getting the following image:

From here I would like to detect the trees.

What I'm doing right now:
First, I'm doing an erosion:

kernel = np.ones(shape=(5, 5),
                 dtype=np.uint8)

image= cv2.erode(image,
                 kernel=kernel,
                 iterations=14)

and getting some overlapping, e.g.:

In order to separate those overlapping contours, I do dilation this time:

image= cv2.dilate(image,
                  kernel=kernel,
                  iterations=12)

Ups, some contours became so small and almost disappeared.
Reversing damage:

image= cv2.erode(image,
                 kernel=kernel,
                 iterations=7)

Ok! Better! But.. all this process seems like a dirty workaround.. Moreover, this won't work for some other images so well (ignore the noise on the left, I should add cropping):

Questions:
 - How do I make sure not to overlap the contours of the trees and not to remove them while dilating/eroding?
 - Is there a better alternative to what I am doing right now?

Edit:
Probably for the full picture I should add that after obtaining a good image with separate contours, I use:

disparity, contours, hierarchy = cv2.findContours(image,
                                                  mode=cv2.RETR_LIST,
                                                  method=cv2.CHAIN_APPROX_SIMPLE)

and plot the contours on the original image:

Edit Nº2:
Adding cross-section of disparity map. Raw data available from here

Plot:

Update on results:

I tried a median filter as it was suggested by @A_A, and got a very good separation of trees from the ground:

ground = medfilt(disparity, kernel_size=125)
trees = disparity - ground
boolean_mask = (trees > 2.) & (trees < 15.)
plt.imshow(boolean_mask, cmap='gray')

The problem with this filter is the speed.
scipy.signal.medfilt took me more than 20 mins to calculate everything. But scipy has also a scipy.ndimage.filters.median_filter. I don't know what's the difference between them but results that I got from it were almost the same apart from some areas close to image borders, and it took somewhere from 2 to 5 mins.


texts: I can't wrap my mind around the concept of filtering and the unit circle ? any help?

texts: Can someone elaborate this theory where there is a link between the stability of a filter and the zeros or nulls being on the unit circle?


texts: Why autocorrelation can be more efficiently calculated using the fft

texts: Can anyone explain why autocorrelation can be more efficiently calculated using the fft ?


texts:
I'd like to calculate the Noise Reduction Ratio of a system.

The system I have is $$ H(z)= \frac{z^{-5}}{1-0.1z^-1}.$$
For FIR Filters this can easily done by $$ \sum h[n].$$

For an IIR Filter this gets more complicated. I thought about making it a recursive formula for h(n) and then calculate some sort of sum which might converge (without any result yet).

Another way would be integrating $ |H(w)|$ by replacing z with $e^{jw}$. But there has to be an easier way. Some sort of trick.

I can imagine two things:

The system can be splited into two or more known systems. The NRR can than somehow be added togehter. (How would one do that?)
$z^{-5}$ is just a delay which shouldn't contribute to the NRR. The remaining part has maybe a well known NRR.

EDIT: UPDATE. I made a mistake

NRR for an FIR Filter is given by: $$NRR = \frac{\sigma^2_{Yv}}{\sigma^2_{Xv}} = \frac{1}{2\pi}\int_{-\pi}^{\pi}|H(\omega)|^2d\omega = \sum_n h(n)^2$$

EDIT 2: What's an NRR

Since I left my textbook in the office I can only show you some definitions from powerpoint slides:

  The noise reduction ratio NRR measures, how much the noise power at the output is less than the noise power at the input
  –
  The last equal sign follows from Parseval’stheorem
  –
  Note that small NRR means large noise reduction
  –
  Small NRR are good NRR

and

  $$\frac{SNR_{out}}{SNR_{in}} = \frac{E[y_s(n)^2]}{E[y_v(n)^2]} \cdot \frac{E[v(n)^2]}{E[s(n)^2]} = \frac{1}{NRR} \cdot \frac{E[y_s(n)^2]}{E[s(n)^2]}$$
  If $x_S(n)$ isn't changed by Filter $H(\omega)$ then:
  $$ \frac{SNR_{out}}{SNR_{in}} = \frac{1}{NRR} $$


texts: How do you determine the power required for a phase shift keying modulation scheme?

texts:
This wiki article has a diagram :

which shows how the message "11000110" is conveyed as a waveform.

I'm trying to understand how the power to transmit this message can be determined.

Do I need to integrate the magnitude of the I(t) and Q(t) signals to determine power? Is there a formula somewhere?


texts:
I cannot tell how decrease/increase in Tx output power should impact EVM?
What is the mathematical relation of I,0 values on the constellation map and output power?

in the context of 802.11, I am trying to understand what is the limiting factor of the output power of a transmitter? (besides power consumption)

The Protocol does not specify output power number in dBm. It specifies EVM threshold for the different rates.


texts: Spectrogram with square or non-square magnitude of STFT: power vs. magnitude

texts:
As seen in this question and answer, to do a spectrogram, it's common to plot either:

the square magnitude $|\text{STFT}(\text{frame}, \text{bin})|^2$ ("power spectrum")
the magnitude $|\text{STFT}(\text{frame}, \text{bin})|$

In an audio context, when to choose the one or the other? For which reason is one more useful than the other?

In the context of noise reduction, there are both Magnitude Spectrum Subtraction or Power Spectrum Subtraction, as seen around page 5-8 in http://dsp-book.narod.ru/304.pdf. When to use the one or the other?

Also, what do most audio editors display? spectrum or power spectrum?


texts: Position Kalman Filter fails to track a constant-acceleration path

texts:
I'm trying to build a discrete Kalman Filter that fuses accelerometer (acceleration) and GPS (position, velocity) measurements. However, I'm finding that my filter can't properly track a constant-acceleration path:

As seen in the top right graph, the difference between the filtered and actual x-positions has a ramping bias, and the variance of the residual noise also increases over time.

My implementation uses the following standard signal model:

$A = \begin{pmatrix}
1 & \Delta t \\
0 & 1
\end{pmatrix},
\quad
B = \begin{pmatrix}
\frac{1}{2}\Delta t^{2}\\
\Delta t
\end{pmatrix},
\quad
C = \begin{pmatrix}
1 & 0\\
0 & 1
\end{pmatrix}
\quad
D = \begin{pmatrix}
0\\0
\end{pmatrix}$

$x_{n+1}=Ax_{n}+Bu_{n}$

$y_{n} = Cx_{n}+Du_{n}$

Where $x = (r_x,v_x)^{T}$, the measurement $y=x$ uses the GPS position/velocity readings, $\Delta t$ is the sample time and the accelerometer reading is used as the control input $u$. My $Q$ and $R$ matrices are both constant, reflecting each sensor's minimum rated accuracy.

Is this poor tracking behavior to be expected, or is there something wrong with my implementation? I also understand that several other alternative signal models can be used for this problem - for example, using the acceleration as a observable state instead of a control input. I've tried implementing this model too, but the filtered result is even noisier than what's shown above. Is one of these models more suitable over the other, and why?


texts:
Now I am trying to design FIR based adaptive filter for rejection of Jamming in GPS device. (just self-learning purpose, Jamming is simple tone)

I've designed NOT-BAD performance adaptive filter for low sampling frequency. (~44100 Hz)

But when I use this adaptive filter for high sampling frequency (16 MHz), transient width of filter is greatly increased.
(I think it is sure because transient width is proportional to 1/degree)

Since it is self learning purpose, I just want to know how big sampling frequency of normal GPS signal for digital filtering is.

If it is not so high as above, I'll quit this design and concentrate on other issues(quantization error etc).

If sampling frequency of normal GPS signal is high as above, I'll study more about that issue.

So, my question is :

How high sampling frequency of normal GPS is ?

Thank you.


texts: PHAT filtering for cross-correlation yields wrong results while simple cross-correlation does OK (Python)

texts:
I work with two signals, the second being perfectly delayed and attenuated, with no reverberation. I implemented cross-correlation algorithm, which yields me these results:

Here is my code for these results:

res1_fft = np.fft.fft(res1_red)
res2_fft = np.fft.fft(res2_red)
prod = np.conjugate(res1_fft) * res2_fft
CC = np.fft.ifft(prod)

abs = [i for i in range(-int(len(CC)/2),int(len(CC)/2))]
plt.plot(abs,np.fft.fftshift(np.real(CC)))
plt.grid()
plt.show()

Please note that I both tried with my data being zero-padded or not. In both case I get a pike at the right delay of 137 samples I imposed.

However as soon as I add the PHAT filtering with this code:

res1_fft = np.fft.fft(res1_red)
res2_fft = np.fft.fft(res2_red)
prod = np.conjugate(res1_fft) * res2_fft
CC = np.fft.ifft(prod)
sign = 1/(np.abs(prod)) * prod
GCC_vrai = np.fft.ifft(sign)

and then plot GCC_vrai, I get weird results, with a pike always in 0:

I also tried this piece of code I found there (with a similar problem seen there too), with a very similar result.

N = len(res1_red)
Ncorr = 2*N-1;  # N size of the signal
NFFT = 2**nextpow2(Ncorr)
R12 = np.fft.fft(res2_red,NFFT)*np.conjugate(np.fft.fft(res1_red,NFFT))
r12_temp = np.fft.fftshift(np.fft.ifft(np.exp(1j*np.angle(R12))))
r12 = r12_temp[int(NFFT/2+1-(N-1)/2):int(NFFT/2+1+(N-1)/2)];

abs = [i for i in range(-int(len(r12)/2),int(len(r12)/2))]
plt.plot(abs,np.real(r12))
plt.grid()
plt.show()

What am I doing wrong?


texts: What are the problems with designing an FIR filter using FFT?

texts:
I'm trying to get an understanding of the relationship between an FIR filter designed from "first principles" using a filter kernel with convolution, and a filter designed in one of two ways using FFT (see below).

As far as I understand, the impulse response of an FIR filter is the same thing as the filter's convolution kernel. (Correct me if I'm wrong.)

Also, in my understanding, the component frequencies (i.e. Fourier transform) of an FIR filter's impulse response is the same thing as the filter's frequency response. And, therefore, the inverse fourier transform will give me back the impulse response.(Again, correct me if I'm wrong).

This leads me to two conclusions (ignoring phase response, or assuming linear phase response):

I should be able to design an FIR filter of arbitrary frequency response by "drawing" my desired frequency response, taking an IFFT to get the impulse response, and using that as my convolution kernel.
Alternatively, I should be able to create a filter by taking the FFT of the input signal, multiplying by my desired arbitrary frequency response in the frequency domain, and taking an IFFT of the result to produce the output signal.

Intuitively, it feels like 1 & 2 are equivalent, but I'm not sure if I could prove that.

It seems like people (and DSP literature) go to great lengths to design FIR kernels with predefined responses, using complicated (to me) algorithms like Chebyshev or Remez (I'm throwing out some names I've read, without really understanding them).

Why go to these lengths, when an FFT/IFFT transform exists for every possible FIR kernel?
Why not simply draw the exact frequency response you desire, take an IFFT, and there's your FIR kernel (method 1 above)?


texts: I have a scenario where I have two OQPSK modulated signals and I need to increase the power of one signal by 3 dBm while keeping the other signal at the same power level. Could someone guide me on how to achieve this in Matlab.


texts:
I am wondering how to shift the center frequency of the below complex data from 0, to say 20 MHz (10 MHz to 30 MHz- 20 MHz bandwidth).  I tried to change the indexing of the bins but I could not get it to shift.  I understand fftshift shifts frequencies about 0 Hz, any way to shift the whole spectrum to something other than 0?

nFFTSize = 64;
subcarrierIndex = [-26:-1 1:26];
nBit = 1000;
ip = rand(1,nBit) > 0.5;
nBitPerSymbol = 52;
nSymbol = ceil(nBit/nBitPerSymbol);

ipMod = 2*ip - 1;
ipMod = [ipMod zeros(1,nBitPerSymbol*nSymbol-nBit)];
ipMod = reshape(ipMod,nSymbol,nBitPerSymbol);

st = [];

for ii = 1:nSymbol

inputiFFT = zeros(1,nFFTSize);

% assigning bits a1 to a52 to subcarriers [-26 to -1, 1 to 26]
inputiFFT(subcarrierIndex+nFFTSize/2+1) = ipMod(ii,:);

%  shift subcarriers at indices [-26 to -1] to fft input indices [38 to 63]
inputiFFT = fftshift(inputiFFT);

outputiFFT = ifft(inputiFFT,nFFTSize);

% adding cyclic prefix of 16 samples
outputiFFT_with_CP = [outputiFFT(49:64) outputiFFT];

st = [st outputiFFT_with_CP];

end

st2 = st*exp(2*pi*20*i*.0000032);

scatterplot(st2), grid;


texts: Amplitude in frequency domain of the same signal measured at different bandwidth (f_max - f_min) and different FFT lines

texts:
Fourier series in complex form

Fourier series in trigonometry form

Here, I'm measuring the same signal (an impulse force - or at least as close as possible to an impulse, represented by a knocking force). The work I'm doing is parameter study (i.e.: studying the effect of different bandwidths and numbers of FFT lines).
This is on the same number of FFT lines but different bandwidth

And this is the same bandwidth and different FFT lines

The main conclusion I can draw is that changing one parameter (while keeping the other is the same) would change the amplitude of the frequency spectrum. I can explain this via math above (smaller frequency resolution leads to longer measuring time, and thus lower Fourier coefficients).
But what does it mean in a practical and realistic manner? As in for actual sampling process, how would changing the parameters (measuring time, sampling points...) affect the frequency spectrum of the same signal?
Edit: The smaller frequency resolution would eventually lead to lower amplitude in frequency spectrum is considered as analytical math by my instructor, while the actual system works with numerical math. So a more "physical" (or "what actually happens") answer is needed.


texts: Intuitively, why is windowing function a low-pass filter?

texts:
I'm trying to intuitively understand why a windowing function of a signal is a low-pass filter. I know its DTFT is a sinc function, making it mathematically a low-pass filter, but my logic says truncating a signal with a windowing function would preserve high frequencies better than low frequencies, since high frequencies have a better chance at showing up in the window.

Thanks.


texts:
I am studying about the KLT tracker algorithm. In the KLT tracker, we can try to estimate the optical flow for corner points (because these are good features to track) between frames and the computed optical flow allows us to obtain the trajectory of the moving object.
This problem is formulated as an LK-alignment problem where we are given an image $I$ and a template $T$. There exist some warping function $W(x;p)$ that maps the image $I$ to the template $T$. The optimisation for LK-alignment is given as $$min_p \sum_x [I(W(x;p)) - T(x)]^2$$ However, this optimisation is expressed as $$min_{\Delta p}\sum_x[I(W(x;p+\Delta p)) -T(x)]^2$$ assuming that we have an initial estimate of $p$.
My question is why can the tracking problem be framed as an LK-alignment problem. I do not see how solving for $p$ or $\Delta p$ gives us the optical flow vectors for each frame.
Would appreciate much enlightenment over this topic.


texts: What level of correlation coefficient suggests correlation?

texts:
In Python (and Numpy), I recently did:-

a = np.random.randint(0,256,1000)

b = np.random.randint(0,256,1000)

and when I calculate the correlation coefficient $r$ using:-

np.corrcoef( a,b )

I get a value of 0.07944437.  Rounded, that's 0.08, 8% or about $1 \over 13$ th.

$r$ is said to typically range -1.0 thru 0.0 to +1.0.  I've never had a 0.0 and don't think that it's possible as random correlations can always be found between any two series. It's like asking "How small is small?" 0.08 doesn't seem that small. Yet clearly $a$ is not correlated to $b$ as both originate from the good quality Mersenne Twister deep inside Python.

So if $r = 0$ is empirically impossible, when does it actually become significant and suggest a real correlation?

Note 1.  There is of course the 68–95–99.7 rule. This covers a $3 \delta$ range of certainty. The atomic guys use $5 \delta$ to confirm a discovery. Do  signal guys do anything similar, and if so, how would they apply it?

Note 2. With the greatest of respect, How do I know quantitatively if the correlation of two time series is significant? seems ambiguous.


texts:
I'm detecting vehicle honking from environmental sounds (engine noise, speech music, siren etc.) based on binary classification. Common audio features (spectral flux, centroid, mfcc, harmonicity etc.) seemed not doing well --- we get model with nice sensitivity but the problem is that the false-alarm rate is high. Failure cases includes sirens and speech.

Any features that describes most common vehicle honking sounds with discrimination from similar sounds for example, siren or speech? Thanks.


texts: Thoretically, what makes multi-pitch detection so difficult?

texts: For multi-pitch detection, I mean for example pressing several piano keys at the same time. For signal pitch detection (pressing one key at a time), after FFT, mostly the frequency with the largest magnitude is desired frequency.  However, for multi-pitch detection, let's just say press two piano keys, the two frequencies with the first two largest magnitudes are not the frequencies of keys I pressed.


texts: Time-Frequency Analysis of Big Data - Data Size Reduction: averaging the most appropriate method?

texts:
Explanation:

I would like to analyse the data from an experiment, which investigates the performance of a mechanical component using sensors, that has generated 2000 CSV files. Each file contains 513 Rows x 1220411 Cols, and they are in spectrogram format (columns are time and rows are frequency):

| Time (s)| 0.0000 |0.000164|...
|:--------|--------|:-------|:-------
| 1.52kHz |  2747  |  350   |...
| 3.05kHz |  2996  |  420   |...
| 4.57kHz |  4078  |  300   |...
|   ...   |  ...   |  ...   |...

I have plotted my 3D chart of the first 100 rows using persp3D():

persp3D(x,y,z, theta=45, phi=5, xlab="Frequency (kHz)", ylab="Time (s)", axes=TRUE, expand=0.5, shade=0.2)

I would like to extract 1-4 column (1 sample of time for all frequencies) from each file to get a table of data with the total of 2000-8000 columns (for 2000 files with all the frequencies) and plot this to get the 3D plot of the experiment.

Question: best data reducing method

I would like to know what is the best method to make sure that those 1-4 columns represent the total data set in each file?
Is averaging simply a good method in this case? What are the alternatives?


texts:
am trying to implement OQPSK modulation with half sine pulse shaping. If I simulate the OQPSK modulation without pulse shaping, the simulated curve is similar to the theoretical value. However, when I introduce pulse shaping, the simulated values are lesser than the theoretical case. I am not sure where I am going wrong. Below is the code where I am using pulse shaping:

  Rc = 2e6;       % Chip rate
  Tc = 1/Rc;      % inverse of chip rate
  Fs = 12e6;      % sampling frequency
  dt = 1/Fs;
  over = Fs/Rc;   % sampling factor
  TSim = 0.01;
  N = floor(TSim/4e-6);
  data = randi([0 1],N,1);

  hMod = comm.OQPSKModulator('BitInput',true);
  modulated_data = step(hMod, data);

   %% Upsample
  over_data = upsample(modulated_data,over);

  %% Pulse Shape
  time = 0:dt:2*Tc;
  hSineflt = sin(pi*time/(2*Tc));
  tx=conv(over_data,hSineflt,'same');

  %% AWGN channel
  hAWGN = comm.AWGNChannel('EbNo',SNR);
  noisy_signal = step(hAWGN, tx);

  %% Matched filter
  rx=conv(noisy_signal,hSineflt,'same');

  %% Downsample
  downsampled_signal = rx(1:over:end);

  %% OQPSK Demodulation
  hDemod = comm.OQPSKDemodulator('BitOutput',true);
  demod_data = step(hDemod,downsampled_signal);

  %% BER
  hError = comm.ErrorRate('ReceiveDelay',2);
  errors = step(hError,data, demod_data);
  BER_sim = errors(1);

This is the part where I calculated the theoretical value

  BER_Theo = berawgn(SNR,'oqpsk','nondiff');

I expect both BER_sim and BER_Theo to be almost the same. But I am unable to figure out where I am going wrong.


texts: How to estimate the frequency characteristics of a nonlinear chirp?

texts:
Given a partial time-domain recording of a noise-corrupted chirp that asymptotically approaches a frequency, and the function that determines the non-linearity, how can I accurately estimate the final frequency it approaches?

Here's an example of what I'm talking about. A radio receiver hears noise until around sample x=5100, then the chirp starts at f=0, and quickly sweeps up to an unknown frequency, but then fades out around x=5500 to some annoying filter ringing caused by a hardware front-end band-pass filter.

X: sample count, Y: ADC count

Ideal solution will have the following characteristics:

Accurately predict the final frequency
Work over a range of final frequencies
Be invariant to different starting frequencies
Be invariant to different stopping frequencies
Relatively quick to compute

Approaches Tried So Far:

Estimating instantaneous frequency using zero-crossing, then curve-fitting. It seems that noise and filter-ringing causes errors in the frequency calculation that make it very sensitive to slight changes in waveform shape. It works better but not well if the signal is resampled to a higher sample-rate first.
Estimating instantaneous frequency using autocorrelation, then curve-fitting using the strongest lags from autocorrelation Function. This works well given knowledge of the incoming frequency so that the autocorrelation window size can be picked ahead of time, but fails when the ACF window isn't wide enough to see the whole signal. If the ACF is too wide, it'll sometimes see a harmonic as the strongest lag.
Curve-fitting the raw signal. This tends not to converge if there's even a little error in the initial conditions, or noise in the signal.
FFT/STFT, then curve-fitting the bins. What I ran into here was that it seemed the bins were destroying the shape of the signal. I couldn't find a good size of bin that gave enough frequency-resolution and time-resolution.

An example recording is available here.


texts: Advantage of using FIR Polyphase Decimator over just computing FIR output every D samples

texts:
Let's say you want to implement an FIR filter $h[k]$ with $L$ taps and downsample the output:

$v[k] = x[k] * h[k]$

$y[m] = v[kD]$

The "naive" way would be to compute all the samples of $v[k]$, then compute $y[m]$.  This should have complexity of $O(LD)$ to process a block of $D$ samples.

The "less naive" way would be to just compute every $D^{th}$ sample of $y[m]$ directly:

$y[m] = v[mD] = \sum_{l = 0}^{L-1} x[mD - l] h[l]$

This should have a complexity of $O(LD/D) = O(L)$ to process a block of $D$ samples.

The "smart" way, according to everything I read (e.g. Proakis & Manolakis Digital Signal Processing: Principles, Algorithms, and Applications, 4th ed. p.771, also here), is to use a polyphase filter/decimator.

This seems to have the same computational burden: every $D^{th}$ timestep, you are computing the output of $D$ filters, each of which has $L/D$ taps, for a total complexity of roughly $O(DL/D) = O(L)$.

Is my analysis correct? If so, what is the advantage of using a polyphase decimator if the "skip $D$" method has the same complexity and is much easier to implement?


texts:
For given 4 conditions, I have to find out what is $x(t)$ with period of 3, and I don't know if $x(t)$ is real or not.

For fourier coefficients $x_k$,
$$1.\ x_k=x_{k+2}$$
$$2.\ x_k=x_{-k}$$
$$3.\ \int_{-0.5}^{0.5}x(t)\, \mathrm dx\ = 1$$
$$4.\ \int_{1}^{2}x(t)\, \mathrm dx\ = 2$$
from condition 1 and 2, I've got this equations

fourier series expansion of $x(t)$: $$x(t)=\sum_{k=-\infty}^{\infty}\ x_k \exp(j{2\pi \over 3}kt)$$
by condition 2,
$$x(t)=x_0+\sum_{k=1}^{\infty}x_k(\exp(-j{2\pi \over 3}kt)+\exp(j{2\pi \over 3}kt))$$
and as $\exp(-j2\pi ft)+\exp(j2\pi ft)=2\cos(2\pi ft)$,
$$x(t)=x_0+2\sum_{k=1}^{\infty}x_k\cos({2\pi \over 3}kt)$$
by condition 1,
$$x(t)=x_0+2\sum_{m=1}^{\infty}\ (x_1\cos(2\pi{(2m-1) \over 3}t)+x_0\cos(2\pi{2m \over 3}t))$$
For condition 3,4 when I integrate $x(t)$, it gives me too complex sums.
$$\int_{-1/2}^{1/2}\cos(2\pi{2m-1 \over 3}t)\, \mathrm dt\ = {3 \over 2\pi (2m-1)}2\sin(\pi{2m-1 \over 3})$$
$$2\sin(\pi{2m-1 \over 3})={\sqrt3, 0, -\sqrt3, \sqrt3, 0, -\sqrt3, ...}$$
so, $\int_{-1/2}^{1/2}\sum_{m=1}^{\infty}\cos(2\pi{2m-1 \over 3}t) \mathrm dt\ ={3 \over 2\pi}\sum_{n=1}^{\infty}({\sqrt3 \over6n-5}+{0 \over6n-3}+{-\sqrt3 \over6n-1})=0.75$
This kind of summation comes out over and over. Of course, calculator can solve it because it converges. But this was problem from my mid-term exam and calculator was not allowed. This thought process(?) complicates the problem too much, so I wondered and asked here if there could be other way to solve this problem much simpler than this.


texts: Selectively improve the frequency resolution of FFT (DFT)

texts:
I have a digitally sampled signal(sampling frequency = 4000 Hz) which contains analog frequencies from 0 to 1000 Hz.
I wish to have a frequency resolution of 100 Hz(like 0,100,200,300 Hz .. 1000 Hz) and I have to improve the frequency resolution at the 0 to 100 Hz range by 10 Hz. That is at the frequencies 0,10,20,30,... 100.
So the frequency components I wish to detect are [0,10,20,...,100,200,300,...,1000].
In other words, I want more frequency resolution at the 0-100Hz range and less resolution at 100-1000Hz range.
How this is possible through FFT or related methods.


texts: python: how to compute the gray level histogram features as mentioned in the paper, and

texts:
Hi I am extracting the grey level features of image mentioned in this paper (part 4.5 low level features).

  "We describe two low level features that are particu- larly important for photo quality assessment – contrast and brightness. Professional photos usually have higher contrast than snapshots. We measure the contrast of an image as fol- lows. First, we compute the gray level histogram Hr , Hg , and Hb for each of the red, green, and blue channels, re- spectively. Then, we compute the combined histogram H, where..."

I can't clearly understand the above part in the paper.
As my understanding about grey level features, take the above image as example.
Use the opencv convert the color image to grey image. as the following code:

image = cv2.imread('lenna.png')
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

and then using the following code to get the histogram with the following code:

  hist = cv2.calcHist([gray], [0], None, [256], [0,255])

I want to ask, and I don't know the meaning of gray level histogram for each channel.

How to we compute the gray level histogram Hr , Hg , and Hb for each of the red, green, and blue channels?
How to compute the the contrast quality, qct mentioned in the paper?


texts:
Assume i have received a known analog OFDM symbol (such the Short-Training-Field (STF) in WiFi systems), call it $s(t)$ with a known DFT per subcarrier, in the presence of multipath and absence of noise.

$$r(t) = \sum\limits_{i=1}^p \gamma_i s(t-\tau_i)$$
where $\gamma_i$ and $\tau_i$ are complex gains and the delays of the multipath, respectively. Without loss of generality, let $\tau_1 > \ldots > \tau_p$.

My question is simply, what is $\tau_1$? I could imagine what $\tau_2 \ldots \tau_p$ are. They are just relative delays with respect to $\tau_1$. But i don't seem to quite grasp how to quantify or understand what $\tau_1$ means. Is it absolute ? Is it with respect to something ?

Any suggestion is much appreciated.

Thanks


texts: Why MFCC have varying length dimensions?

texts: I have recently started to use MFCC as a feature for a project. However, I am unable to understand why MFCC have varying 2nd dimension size even when I chose the number of coefficients to be 40.


texts:
I am generating a time-stream data set that needs to be filtered.  The data-set is pretty large (i.e. too large to be filtered in one go).  I'm trying to split the filtering operation up into chunks but I obviously get a spike at each chunk start.  Is there a correct way to split a large data set up into smaller ones and have a continuous result - i.e (pseudcode).

timestream # variable with say 1e9 data points
b,a = sig.butter(1, 0.16, 'lowpass')

Split the timestream up in to say 1000 "chunks"

f_t = []
for c in chunks:
   t_data = sig.lfilter(b,a,c)
   f_t.append(t_data)


texts:
A signal sampled in time domain has replicas of its spectrum in frequency domain.
This implies that to transmit a sampled signal we need ifinite bandwidth.How is this possible?


texts: Why is my discretized transfer function unstable when my discretized state-space model is stable?

texts:
I am working with a 30-state, 14-input linear model that is described by a state-space model:

model_state_space = ss(A, B, C, D);

The model is extremely slow (it describes the thermal behavior of a building), so when I discretize it I choose a sample time of 900 seconds and plot the step response from one of the inputs:

[y,t] = step(model_state_space);
plot(t, y(:,:,end))

I obtain a similar plot when I compute the step response of the discretized state-space model, or of the equivalent transfer function continuous model:

hold on
[y,t] = step(c2d(model_state_space, 900));
plot(t,y(:,:,end), 'r')
[y,t] = step(tf(model_state_space));
plot(t,y(:,:,end), 'k')

But when I compute the impulse response of the discretized transfer function, the system appears unstable:

[y,t] = step(c2d(tf(model_state_space), 900));
plot(t,y(:,:,end), 'g')

Can someone please explain why this is so?


texts: Time domain representation of a modulated signal when symbol rate >> carrier frequency?

texts: I need to understand the relationship between carrier vs. bits with PSK modulation, and what the signal in the time domain would look like if my symbol rate is far greater than my carrier frequency. How can I put more symbols on a carrier which is at a set frequency? For example how would I put 1000000000000000 bits on a 1 Hz carrier wave, is this possible and what would that look like?


texts:
I'm trying to do shot boundary detection using SVM's.

I'm differentiating between 3 shots, cut, fade and dissolve + (normal frames).

I had used two features, shannon's entropy and variance, it gave me 40/635 boundary frames correctly predicted. That is poor accuracy.

What features can I use to improve this?

Training videos: 2 videos from TRECVID (~50000 frames)
Test Videos: 1 Video TRECVID (~13000 frames)

EDIT:

So one must be careful in choosing features.

So I fixed a problem that was occuring, if anyone comes across this post, then the above features are okay and will give you pretty good results for CUT detection.


texts:
Let us assume that I wouldd like to create a digital lowpass with 5 taps. I am not using a window function but just the impulse response from my Fourier integral.

I get the following coefficients: $h=[0.5\, 0.86 \,1 \,0.86\, 0.5]$.

What I don't understand: how does the filtering work? Assuming that I have a constant signal  $x$ with '1' and I want to use the filter on that.

I applied the filter and what I get is

$$y=[0.5 \,
1.36\,
2.36\,
3.22\,
3.72\,
3.22\,
2.36\,
1.36\,
0.5
]\,.$$
But now i am clueless how to draw this? The output is bigger then my input signal?


texts:
I am trying to perform motion compensation for an airborne stripmap SAR data following the book Synthetic Aperture Radar Processing by Franceschetti and Lanari. The algorithm flowchart is on page 145. My questions are:

Can I use the same phase equation given for 1st order MoCo for range resampling?
What is the right domain for performing these compensations? My understanding is that the range resampling can be done in frequency domain to achieve a shift in time domain and the phase compensations can to be done (via multiplications) in time domain for removing the phase error $ e^{j \frac{4 \pi \delta(r)}{\lambda}}$.


texts: Is there a way to double all frequencies?

texts:
Let's imagine I have an input signal that looks like this:

$$
f(t)=\sin\Big(2\pi(2)t\Big)+\sin\Big(2\pi(3)t\Big)+\sin\Big(2\pi(5)t\Big)
$$

And I want an output signal that scales each frequency proportionally by some factor. Here, I scale everything by a factor of '2':

Is there a way to do this?  The frequency domains look like this:


texts: Forming the estimation manifold for a 2-D beamformer

texts:
Given an array of L sensors in a uniform planar array, with known positions in space, I would like to estimate both the angles of arrival and delays from a point source under multipath conditions simultaneously, initially using a Bartlett beamformer. The data from the array is in the frequency domain, taking the form of an L x Nf matrix, with Nf bins around a carrier frequency f.
Based on this paper, I have an L x M angular response vector a, and an Nf x K delay response vector d, where M and K are the resolutions in angle and delay respectively. From these, I wish to create an estimation manifold (defined as U in the paper), which will become the 2-D analogue of the array response vector in the case of 1-D beamforming.
From equation 1 in the paper, it appears as if I need to set M = K, with their product becoming the estimation manifold U, however this produces a level surface when plotted. Alternatively, letting M != K and taking the Kronecker product means this manifold will quickly grow beyong a reasonably storable size, especially when M, K or L are large.
I have also tried slicing the input data into L x 1 slices and treating this as a 1-D estimation problem, which was then projected into 2-D using the delay response vector d, however as this removes the temporal relationships, this did not yield relevant results either.
Is there a better way of forming this estimation manifold, and eventually applying it to a Bartlett beamformer?
Also, for full disclosure this is part of university coursework.


texts:
as the definition of normalized frequency is the frequency divided by the sample frequency, I'm confused about my signal so as an example a signal sampled at ($f_s=500 000$) Hz when applying a FFT, an impulse at frequency of $87$ Hz occurred(true value confirmed) while it occurred at $0.117$ Hz in normalized frequency so we know that 1 for normalized frequency is sampling frequency($f_s$) but $0.117$ Hz normalized frequency how it calculated according to $f_s$ to get $87$ Hz while it's not $87/(f_s/2)$!!
I will be grateful for your help


texts: Cut-off frequency in Fourier-Spectrum of a noise message

texts:
I have a problem involving cut-off frequency.

I have a voice message in an audio file .mat which have an eletromagnetic noise. I am trying to recover the original message by applying digital Butterworth filter with MatLab. To this, besides of design the Butterworth filter, also I have to find the cut-off frequency to eliminate all, or partially, the electromagnetic noise.

The picture below shows the Fourier Spectrum of the voice message with the electromagnetic noise. How could I find the cut-off frequency to eliminate the noise? Thanks in advance!

Sampling frequency = 44.8 kHz.


texts: May I know how to find Time frequency power spectrogram of a signal?

texts:
I am analyzing EEG signal and I am unable to understand the term Time-Frequency Power Spectrogram. Please Can anyone help me?

Thank you in advance.


texts: What causes High Intensity Noise after Acoustic Echo Cancellation for a VoIP?

texts:
I am trying to implement a Acoustic Echo cancellation algorithm in Java for a VoIP system. This AEC will help to cancel the echo (coming from Speaker) in microphone input at the same PC. The sound frequency is 16 kHz mono and 16 bits PCM frame. At first I tested my codes with artificial sound stream. I use a sound file as speaker output, and then put 0.6 second delay on the same sound file to simulate microphone input. My AEC works well with and without double-talk.
My problem is: my AEC can not work when I test it in a real scenario including two remote PCs. At the first a few minutes, there are a lot of high intensities noise. After I transfer the 16-bit PCM frames into SHORT in Java, these noises almost reach the maximum value of SHORT type.  These noise will go down after a few minutes, but echos are not cancelled completely. I have tried to adjust the size of FFT (currently, it is 256),  and the length of filer tail (Currently, it is 9600 frames for a 0.6 second standard echo delay in 16 kHz ) for a lot of different configurations, but none of them works to me.

Does someone have any clue or ideas on this problem? My implementation is mainly based on the two papers below.

J. S. Soo, K. K. Pang Multidelay block frequency adaptive filter, IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-38, No. 2, February 1990.
Valin, J.-M., On Adjusting the Learning Rate in Frequency Domain Echo Cancellation With Double-Talk. IEEE Transactions on Audio, Speech and Language Processing, Vol. 15, No. 3, pp. 1030-1034, 2007.

Additionally, there is an another question about my artificial testing data. If I use the sound of human being talking to simulate speaker output and microphone input, the 'echo' can be cancelled significantly. But the AEC can not converge completely if I use some repeating simple sound (like, calling sound of Skype). Initially, I thought that the bottom one would lead to a quick convergence, but it is not. Maybe this is relevant to my question about, but I am not sure. Also is there someone having any opinion on this?
(I tried to post a image about these two sound types, but I do not have enough reputation to do that T_T, sorry for this)


texts:
I've run into an issue fiddling with the FFT/IFFT routines given here (http://www.dspguide.com/ch12/3.htm). I've rewritten the code to C++ and played with it. The idea is to read a raw wave file (encoded as 32 bit floats),
find its FFT (complex) and then use that FFT as input to the InverseFFT function.

I expected to get the original sound data (samples) and I do for FFT lengths of up to 2048. For greater lengths (>= 4096 ), the code below doesn't (seem to) work. As you can see, my test consist solely of checking values
at indexes 20,30 and 40 for both the original samples and the result (inverse FFT) array. I've tried this on several files (created with Audacity) with the same results.

    int main(int argc, char *argv[])
{
    unsigned fs;
    float *Ix;
    QCoreApplication a(argc, argv);
    cout << "opening file...";

    char *pData= OpenRaw("C:\\test\\250_350_c.raw",fs);

    float* Rx = (float*)pData;
     cout << "original data { [20,30,40] } " << Rx[20] << " " << Rx[30] << " " << Rx[40] << endl;
    const unsigned numSam = 2048;
    Ix=new float[numSam];
    for(int i=0;i<numSam;i++)
        Ix[i]= 0;

    FFT(Rx,Ix,numSam);

    float *Ifft = InverseFFT(Rx,Ix,numSam);

    cout << "IFFT data { [20,30,40] } " << Ifft[20] << " " << Ifft[30] << " " << Ifft[40];
    return a.exec();
}
char* OpenRaw(const char *fileName, unsigned& fs) {
    char *pData;
    unsigned fileSize = 0;
    int bRead = 0;

    ifstream file (fileName, ios::in | ios::binary);
    if (file.is_open())
    {
        cout << "reading file..." << endl;

      file.seekg (0, ios::end);
      fileSize = file.tellg();
      cout << "file size " << fileSize << endl;

      pData = new char[fileSize];
      file.seekg (0, ios::beg);
      bRead = file.readsome(pData, fileSize);

      if(file.fail())
          cout << " error reading file..." << endl;
      else
          cout <<  "bytes read.." << bRead << endl;
      file.close();

    }else{
        cout << "error reading file";
    }
    fs = fileSize / sizeof(float);
    return pData;

}
float *InverseFFT(float *Rx, float *Ix, int fftSize){
    float *pWaveData = new float[fftSize];
      for(int i=0; i<fftSize; i++){
          Ix[i] = -Ix[i];
      }
      FFT(Rx,Ix,fftSize);//calling the FFT function
        for(int i=0; i<fftSize;i++){
         pWaveData[i] = Rx[i]/fftSize;
         Ix[i] = - Ix[i]/fftSize;
        }
    return pWaveData;
}

I haven't posted the FFT code here because I've used that code in the past and it always worked OK for me (I can still post it here if someone wishes to experiment with the code).

Does anyone know what might be causing this behaviour?

UPDATE 1: I've implemented the second method found in the link given by Peter below, but the result is still the same (the IFFT works for fft lengths of only up to 2048). My original IFFT method (given above) is actually the method #4 found in that same link.

float *InverseFFT(float *Rx, float *Ix, int fftSize){

FFT(Rx,Ix,fftSize);
float *pWaveData = new float[fftSize];

  for(int i=1; i<fftSize/2; i++){
      pWaveData[i] = Rx[fftSize-i]/fftSize;
      pWaveData[fftSize-i] = Rx[i]/fftSize;
  }

return pWaveData;

}

UPDATE 2: After googling some more, I found this ( https://www.mathworks.com/matlabcentral/newsreader/view_thread/247760?requestedDomain=www.mathworks.com ).
Well, unfortunately, it seems that longer fft lengths may adversely influence an FFT/IFFT algorithm, especially if it's a power of 2 based algoritm.

Thank you


texts: How transmission speed and bandwidth are linked?

texts:
I don't understand why if I have a larger bandwidth I can transmit data faster.
Is this linked with this property of fourier transform?

thanks.


texts: Optimal inference for nonlinear state space models

texts:
When considering a linear-Gaussian state space model, it is often referred that, optimal inference is tractable which is very rare in state space models. When considering a nonlinear state space model, it is said that, in general the posterior distribution is intractable, therefore, inference methods should be used to approximate the posterior.

What I would like to ask is that, what is the exceptions for nonlinear state space models? Which nonlinear state space model is tractable analytically? Are there any examples?

PS: For instance, Doucet says:

  Given the ﬁltering distribution one can then routinely proceed to
  ﬁltered point estimates such as the posterior mode or mean of the
  state. This problem is known as the Bayesian ﬁltering problem or the
  optimal ﬁltering problem. [...] Except in a few special cases,
  including linear Gaussian state space models (Kalman ﬁlter) and hidden
  ﬁnite-state space Markov chains, it is impossible to evaluate these
  distributions analytically.

To reexpress, I would like to ask that, which nonlinear state space model is tractable?


texts: Matlab: Correct method of generating Rayleigh fading channel

texts:
A Rayleigh flat fading channel has 1 tap and so there will be no intersymbol interference.
if

N= 1000;
 Then,
h=(1/sqrt(2))*(randn(N,1)+1i*randn(N,1));

will generate 1000 coefficients which are complex Gaussian of zero mean.

Question 1: Is this a multi-path Rayleigh channel of number of taps = 1000?(Example taken from : http://www.mathworks.com/matlabcentral/fileexchange/30477-ber-curve-for-bpsk-in-rayleigh-channel/content/BER_BPSK_Rayleigh.m )

But, there are models such as Clarke's method of sinusoid and Jake's algorithm for modeling the Rayleigh channel as FIR. If the method in Question 1 and these algorithm create multi-path Rayleigh, then when do we use the Algorithm?

Question 2: In the code and explanation therein for Rayleigh flat fading code , a channel of length N is created but it is called flat. I am really confused as to what is the proper way to create multi-path Rayleigh fading and only flat Rayleigh fading (1 tap).

Please help.


texts: GMSK Big error in demodulation MatLAB Simulink

texts:
I am trying to do a PLC simulation system, but my problem is in GMSK demodulator, in my simulation in Simulink MatLAB my output signal is very different from the input signal, below the model:

Below the block parameters:

I start doing this simulation yesterday, the values of each block is the default, and even without de AWGN block the error still appears. Someone knows what is happening?


texts:
I am studying linear minimum mean squared error (MMSE) receivers. The technique used in it suggests to make colored noise white by multiplying it with the invertible covariance matrix $K_z^{-\frac{1}{2}}$, with $K$  given by

$$K_{zk} = N_0 I_nr +\sum_{i≠k}^{n_t}P_ih_ih_i^*$$

where $h_i$ and $h_i^*$ are the channel gain and channel gain conjugate transpose, respectively, and $n_t$ is the number of transmit antennas. The colored noise is because of interference from other users (I guess so!).

I am confused as how the multiplication of $y = hx + z$ with $K_z^{-\frac{1}{2}}$ can transform colored noise to white though I can guess that may be multiplying with inverse cancel out the colored part and the noise is white now. Am I correct?


texts: How to make a definite integral in python?

texts:
Helo!
I have this exercise to do but I don't understand how to define the integral.
$$ v(t) = \int_{-\infty}^{t} y( \lambda) d \lambda $$ considering $$y(t) = 10 \Lambda \left( {t-1 \over 1} \right)$$
Can someone help me please?


texts:
I want to study Detailed explanation of hybrid systems?Which incorporate both continuous and discrete time signals & systems?
for example
In which a continuous-time input signal is transformed into a discrete-time output signal and vice versa. An analog to digital (A=D) converter is an example of a hybrid system.
Where can i find detailed explanation to quench my curiosity/thirst of this topic , i have tried searching google but no useful result


texts:
I'm working on a object tracker based on KLT algorithm using matlab.
Can anybody guide me how I can find the features and calculate the distance in the sequence of the frames?
Thanks.


texts: How to prove that this impulse response is stable?

texts:
This is the impulse response:

Can anyone give a detailed method to prove that it is stable or not ?


texts: Estimating frequency when there's little isolation between signal and noise frequencies

texts:
Say I have a signal which is guaranteed to have a frequency between 110-120 Hz but is corrupted by interference signals that're very close to this frequency range. For example, let the interference signal have a frequency of 150 Hz and an amplitude of 10% if the fundamental.

What kind of filter or filtering technique would be suitable for estimating the frequency of the fundamental while ensuring that the group delay is low? I'm typically looking at a delay between 1-2 cycles of the fundamental frequency. I've experimented with Windowed sinc FIR filters having a cutoff frequency of 125 Hz & a stop band edge of 130 Hz but such a narrow pass band results in a high amount of group delay which'd be unacceptable for this application. The sample rate is pretty high at 25kHz.


texts: What is a good FFT window function to reject DC?

texts:
I'm using an FFT to analyze what is essentially the power envelope of a signal (see here for info on the containing project), and, since power numbers are always positive, to eliminate the DC component I'd like to use a window function that is 50/50 positive and negative, vs the usual all-positive function.

I've taken the "flat top" function, removed the a0 bias and converted it from cosines to sines, but I'm not sure that's optimal (or even meaningful).

Any suggestion?


texts: What exactly is the effect of a Hann Window on the FFT output?

texts:
In the code I am trying to understand, there is an usage of Hann window multiple times (There's overlapping frames, the size of each is 512 and the overlapping factor is 75% so each time we capture  128 samples we build a new frame by shifting the previous one by 128 to the left). For each Frame we apply a Hann window before the FFT, so the first 128 samples were actually "framed" 4 times.

However I don't understand the purpose of this windowing. Can someone please tell me what its effect is on the FFT ?


texts: Length of DFT defines whether I see the harmonic or the dominant frequency

texts:
I have a number of signals that are periodic. I use an fft transformation to obtain the dominant frequency of each signal. In order to increase the frequency resolution I zero pad the signal before the fft. However I noticed that when the frequency is low, the peak of the fft corresponds to the first harmonic. A sample of my signal looks like this:

My code is the following:
N = 4000; %signal length
minf= 0.5; % Hz
maxf= 20; %Hz
resf = 0.01; % Hz resolution
fs = 2000; % Hz sampling frequency

NFFT = fs/resf;
Freq = fft(mysignal,NFFT);
Freq = Freq((minf/resf)+1:(maxf/resf)+1,:);

In this case the fft will perform well and detect the 1.12Hz frequency. However I noticed that when my signal's frequency is below 1Hz, as dominant frequency I will get the 1st Harmonic as the pick of the fourier domain plot and not the dominant frequency. In that case I can change the line  NFFT = fs/reso_freq; to NFFT = fs/reso_freq/2; which will actually give me the correct frequency, but it will not work for when my signal is above 1Hz. Can someone explain to me why this is the case and if there is a universal way of finding the correct frequency without knowing what to expect beforehand.
Thank you


texts: 2D Gaussian Seperation into 1D Gaussian components

texts:
I have been reading this http://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm and I am having trouble to perform the seperability.

So I kinda did it in paper. Lets say y Gaussian function is G(X,Y), then seperating them will become G(X)G(Y), and then I will need to calculate the 1D component for X and 1D component for Y. Then I can pass over my image twice using the two components each time.

I do have a couple of questions though (one of them is more general):

1 - When I calculate the 2D Gaussian why do I need to normalize it by dividing each pixel of the template with the sum of the template?

2 - Do I need to do the same for the 1D components aswell?

3 - In the link it says:

  Figure 4 shows the 1-D x component kernel that would be used to
  produce the full kernel shown in Figure 3 (after scaling by 273,
  rounding and truncating one row of pixels around the boundary because
  they mostly have the value 0. This reduces the 7x7 matrix to the 5x5
  shown above.). The y component is exactly the same but is oriented
  vertically.

Why scale the 1D component by 273? 273 is the sum of the 2D template not the 1D. And how does he truncate 1 row? The 1D X component is 1 row by itself...

Thanks
I hope my questions make sense.


texts: How do i calculate the accuracy of an algoruthm?

texts: I'm doing research on an machine learning image processing algorithm (segmentation algorihtm)i.e. computer vision performs when an image takes as an input and objects in an image are segmented and counted. i wanna ask how can i calculate the accuracy of my model in MATLAB. i searched many links, there's have some training taken to get the accuracy but i'm confused in it. i'm working in machine learning algorithm not deep learning. so is there some other method to find accuracy except it or not?


texts:
I am given the following discrete time transfer function : $$G_d(z^{-1})=z^{-d}\frac{b_0+b_1z^{-1}}{1+a_1z^{-1}}$$ which has the following impulse response $$g_d[n]=\{0,1,-0.1,-0.05,...\}$$

How can I find the transfer function?($b_0$, $b_1$ etc)


texts: Recreating an image from random samples?

texts:
Assume that I have samples from an image, but the samples aren't on a grid, but instead there are multiple lowres versions of the original where the sampling locations are slightly jittered, so the pixels in the lowres images haven't been sampled from the exact same location (i.e. the lowres images aren't identical, however assume I now the precise sampling location for each pixel in the lowres images). What methods exist for reconstructing the original image?

The most general case of this is assume you have random samples from the image (making the assumption here that the image is a continuous signal). Are there optimal algorithms for reconstruction? For example, this assumes that the sampling density can vary spatially in the image, so that some regions have very few samples while others are densely sampled.

My DSP books seem to concentrate on methods where some signal has been sampled at fixed periodic intervals, so they don't seem to be of much use.

EDIT: Pointers to literature that discusses this would be great, even knowing what such reconstruction problems are called!


texts:
Im new to Matlab. I need to split the time dimensions of an audio track into frames of 10mS and for every 10mS frame have to compute the FFT to generate the spectrograph.
I need help on this.


texts: BER of QPSK modulation in Rayleigh fading channel model

texts:
i am having trouble deriving a BER analytical expression for QPSK modulation system operating in a Rayleigh Fading channel model.
Could anyone direct me to a reference that would help in this matter.

Thank you


texts: Mid/Side encoding/decoding in frequency domain instead of time domain

texts:
In the time domain let's say we have a signal $L_t$ and $R_t$. Then we can do mid/side encoding/decoding like so:

$$M_t = \frac{L_t + R_t}{2}$$
$$S_t = \frac{L_t - R_t}{2}$$

$$L_t = M_t + S_t$$
$$R_t = M_t - S_t$$

Let's say the we convert a block of $N$ values from $L_t$ and $R_t$ to the frequency domain by zero padding and taking their FFT to yield $L_f$ and $R_f$. Is it possible to do a similar encoding for mis/side in the frequency domain i.e. convert $L_f$ and $R_f$ to $M_f$ and $S_f$? So that the inverse FFT of $M_f$ and $S_f$ would yield $M_t$ and $S_t$?


texts: What is the set of most predictable or in other words most self-correlated signals?

texts:
Which signals or set of signals is the most predictable or in other words most self-correlated (autocorrelation)?

I think one of them is $x(t)=c$, for all $t$, where $c$ is a scalar constant. The other one can be $x(t)=c$, for $t=t_0$ and $x(t)=0$, otherwise (a discrete signal with only one value). $t_0$ is a constant.

Here the book says: "the most predictable signal" is the sum of discrete cosine functions with random phases. Why? My two examples above seem to me more predictable!


texts:
I'd like to find the "best" fit of an ellipse to contiguous, possibly concave shapes such as:

What have I tried?

I thought that one could assign the direction of the major and minor axes $\vec a, \vec b$ of the ellipse by mapping the pixel values to coordinates, mean subtracting, and saving the largest two eigenvectors from a PCA. The seems to work rather well in finding the direction:

My problem is determining the length of these two vectors. For now, I've used $\sqrt \lambda_1$ $\sqrt \lambda_2$, from the eigenvalues of the PCA. This seems to underestimate the length. How can I determine $|\vec a|$, and $|\vec b|$ or alternatively, best fit an ellipse to these shapes?


texts:
I'm doing a rigid body between point clouds using SVD.
Sometimes the routine produces an incorrect rotation matrix with a Determinant=-1, ie a reflection matrix.

Any idea why ?
Is there a valid way of converting the result into a correct matrix?


texts: How do I resample an image to a rotated grid?

texts:
I have an image, I, sampled on a uniform grid:

$\ x_i = i*\Delta x, y_j = j*\Delta y, $

I need to resample this image to a grid rotated counterclockwise by an angle $\ \theta$ around $\ (x_0,y_0)$:

$\ u_i = i*\Delta u, v_j = j*\Delta v $

$\ x(u_i,v_j) = x_0 + cos(\theta)*u_i - v_j*sin(\theta)$

$\ y(u_i,v_j) = y_0 + sin(\theta)*u_i + v_j*cos(\theta)$

How do I best do this?

I understand that I can e.g. do bilinear interpolation to find I from the 4 (x,y) corners around each $\ (u_i,v_j)$.

However I am not free to choose $\ \Delta u $ and $\ \Delta v $ as I wish.

What is the Nyquist criterion for this rotated grid?

Am I right in assuming that it will be?:

$\ \Delta u = cos(\theta)*\Delta x + sin(\theta)*\Delta y$

$\ \Delta v = -sin(\theta)*\Delta x + cos(\theta)*\Delta y$


texts: Why do we always characterize a LTI system by its impulse response?

texts: Why do we always characterize a LTI system by its impulse response and not by another response, like the step response? What does the impulse response have that is so special?