texts: What are some free alternatives to SIFT/ SURF that can be used in commercial applications?

texts:
As far as I understand, both SURF and SIFT are patent protected.
Are there any alternative methods that can be used in a commercial application freely?

For more info on the patent check out:
http://opencv-users.1802565.n2.nabble.com/SURF-protected-by-patent-td3458734.html


texts: How is loop information encoded in FLAC files?

texts:
I'm trying to convert .wav files (samples for a sampled instrument) into FLAC format.  The .wav files have loop information in the smpl chunk.  I understand that various sample players support loops in FLAC files.  However, I can't find a converter that converts .wav to .flac and preserves the loop information.  I could code it myself (say, by postprocessing results from SoX), but can't find any information on how loops are encoded in FLAC files.

I'm using this for sampled instruments with hundreds of .wav files and so need to be able to do batch processing.

Sorry if this is the wrong forum.  It's not quite stackoverflow or sounddesign.


texts: Find Syncword (pattern) in signal, using Correlation

texts:
I have a signal which always has the same syncword ([93 0b 51 de] in hexa).
I also have a signal with additive noise(attacing it, and I want to find my syncword in that a noisy signal using corelation(actually using xcross corelation), I attach the plot of my signal with additive noise(actually the demodulated signal that I de modulated it), my syncword is hidden here with noised signal and other data(like crc ..but doesn't matter for me and Im just looking for my syncword)

those data Im getting from transmitter that's transmiting on frequency 868MHZ , Im implicitly the receiver and sampling in frequency 868MHz , those IQ samples data has packets and Im seeing them on my matlab

those IQ samples data(sampling them on 20480000) I demodulated them by two matched filters and the data the I got is ranges between 1 , -1 . -1 means 0 in binary and 1 means 1 in binary.. (the plot of my demodulated data see the first photo that I attached above..
code in matlab:
>>fs=204800;
>>dk=loadFile('fm.dat'); %Im loading the IQ samples that I sampled them by my hardware sampler , once I do plot(real(dk)) I see the packets in matlab

>> t20=[-40000:40000]*(1/fs);
>> mf1200 = 100*exp(i*2*pi*t20*100000);
>> mf2200 = 100*exp(i*2*pi*t20*-100000);
>> d12 = conv(dk,mf1200,'same');
>> d22 = conv(dk,mf2200,'same');
>> ddif=((abs(d12)-abs(d22))./(abs(d12)+abs(d22)));
>> xCorrResult=xcorr(ddif',random);

Then I used xcorr, but I was a little confused by that. Can you help me? It looks like that signal was too noisy so that no syncword was found in my demodulated data that ranges between (-1 1) ?

the packet that I transmit is in frequency 868MHZ, 50KHZ symbol rate , deviation 100KHZ , and my syncword/data are attached in the photo below:

the file can be opened on matlab by this matlab code (plot(real(y) will give you plot of packets like nullPacketDatanullPacketDatanullPacketData):
function y = loadFile(filename)
%  y = loadFile(filename)
%

% reads  complex samples from the rtlsdr file
%

fid = fopen(filename,'rb');
y = fread(fid,'uint8=>double');

y = y-127.5;
y = y(1:2:end) + i*y(2:2:end);


texts:
I am trying to implement the symbol timing recovery loop and for that I have chosen two algorithm to compare my results one is feedforward and Gardner TED. One thing that is confusing me is smoothing of error after finding timing error and finding the right sampler.

Should I use the sampler after matched filter or before matched filter?
how the new timing is computed from the timing error?


texts: Why many DACs uses series of half band filters for interpolation instead of single one?

texts: And why, according to some sources, these filters have different number of taps (the last is the shortest one)? Does it really reduce a cost of computation?


texts: why is it so hard to simulate up and down position on headphones?

texts:
i'm wondering how come simulating front,back,left and right source position on headphones is easy but up and down is not ? (maybe even impossible ?)

is it due to the shape of the ear ? but should we not be able to modify the input signal according to the shape of the ear and then simulate ups and downs movements ? why binaural impulses don't simulate correctly ups and downs on headphones ?

thanks

Jeff


texts: Why Is the Result of the Convolution of a Row and a Column a Two Dimensional Matrix?

texts:
I'm trying to understand a book chapter on "Algorithms For Efficient Computation of Convolution" (Also on Doc Droid or Scribd) and I know that when calculating the convolution of an image represented by two dimensional matrix (array) with a mask matrix the result is the treated two dimensional image.

How to calculate the convolution of two signals represented by a row and a column?

The chapter I'm reading says that:

  Separable convolution kernel must fullfil the condition that its matrix has rank equal to one. In other words, all the rows must be linearly dependent. Why? Let us construct such a kernel. Given one row vector

  $$
\vec{u}=\left(u_1, u_2, u_3, \ldots, u_m\right)
$$
  and one column vector
  $$
\vec{v}^T=\left(v_1, v_2, v_3, \ldots, v_n\right)
$$
  let us convolve them together:
  $$
\vec{u}*\vec{v}=\left(u_1, u_2, u_3, \ldots, u_m\right)*\begin{pmatrix}v_1\\v_2\\v_3\\\vdots\\v_n\end{pmatrix}=\begin{pmatrix}u_1v_1&u_2v_1&u_3v_1&\cdots&u_mv_1\\
u_1v_2&u_2v_2&u_3v_2&\cdots&u_mv_2\\
u_1v_3&u_2v_3&u_3v_3&\cdots&u_mv_3\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
u_1v_n&u_2v_n&u_3v_n&\cdots&u_mv_n\end{pmatrix}=A\tag{4}\\
$$


texts: Which one can be accompanied by linear filters?

texts:
I have a matrix :
\begin{bmatrix}
1 & 2 & 3\\
1 & 4 & 5\\
2 & 6 & 7
\end{bmatrix}
After doing operation 1 ,
I get
\begin{bmatrix}
0 & 0 & 0\\
1 & 2 & 3\\
1 & 4 & 5
\end{bmatrix}
However operation 2 on same matrix will give :
\begin{bmatrix}
0 & 1 & 2\\
0 & 3 & 4\\
1 & 5 & 6
\end{bmatrix}
Clearly, operation 2 reduces each matrix value by 1, while operation 1 shifts one row down.
Now I see that operation 1 is simply subtraction and thus can be applied using linear filter. But I am not using to linear filters on matrix. Can anyone guide me through this?


texts: Bandwidth of a modulated signal in analog domain

texts:
If i am using QPSK modulation with data rate 4 Mbps with rolloff factor $0.5$, then the bandwidth will be symbol $\rm rate\times (1 + rolloff\  factor)$. So here it will be $\rm 2\ MHz\times 1.5 = 3\ MHz$
This is baseband data bandwidth.
So should I use analog filter of bandwidth $\rm 3\ MHz$ or should it be more? or should it be more than that?


texts: Can a Goertzel filter response be mapped to a Direct Form II Biquad IIR Filter

texts: I am looking at the structure of the Goertzel Filter, and it seems like I can map its coefficients to the standard implementation of the CMSIS DSP Biquad Cascade IIR Filter. Is this possible or should I give up and implement my own? The benefit here is the CMSIS Library can be had in a certified form for my application (I believe it just got rolled in). I believe I can compose the Goertzel response by cascading the Biquad with a Comb network, is this correct? I suppose my main point of confusion is centered from my inexperience in matters such as these, and if the Transposed Direct Form II of the cmsis biquad can demonstrate total equivalence to the Direct Form II suggested by this particular paper [edit]


texts:
I have a device (Large as typical PC casing) with a motor. There can be a condition that the device motor load is unbalanced and the device starts to move around the table and I need to detect this. For this, I have an accelerometer.
But the problem is that I can't come up with a reliable algorithm to detect these vibrations.
Here are sensor data. The solid line is 5 point moving average. The dotted line is raw data. Time is in ms.

I tried to do some averaging, max and min detection and some other simple algorithms but it didn't work, but there is difference between plots, "Vibration" plot looks more periodic and consist of one frequency.
Can you please suggest some algorithms or ways to deal with such a problem in the time domain?
Maybe there is some neat way to differentiate these two plots?
I am using 8bit MCU to read data and detect vibrations so I can't do real-time frequency domain analysis.
Thank you
PS.
I will add plots I made in the frequency domain, maybe this will help.
I took 5 points moving average. Then used a cosine hamming window and took DFT of 64 samples.
Here are plots
Vibrations

No vibrations


texts: Computer vision, Coupon code detection, sugested algorithms

texts:
The telecommunication operator in my country has coupons that used to recharge the balance

I want to develop an Android application, that with the camera can take the scratched number and recognize and automatic recharge.

I think finding the number box, and latter preform an OCR or a pattern match in because are numbers only

I'm new in computer vision I set the environment and test some opencv ups, but I'm losing witch algorithms or sequence of algorithms may use to achieve that I want Exist segmentation, feature detection, blob detection, ... template match but the image maybe rotates scale because the use of the camera

And I'm lose I wonder if exist one algorithm to find the 4 extreme points of a coupon in one image from here I know how to proceed until the OCR part I don't want the code to do it of course

What I want is a simple clue, answer like you may use X algorithm to have Y result

Thanks and sorry for my english


texts: Spectrum Analyzer - Acquisition time computation

texts:
Using a spectrum analyzer, is there a theorical equation that links acquisition duration and acquisition parameters ( Span, RBW, VBW, sweep time and averages number )
Does this equation depends on the equipment or is a general rule?


texts: Why does a phase shift in the frequency domain not produce a time delay in the time domain?

texts:
Suppose I have an input signal which represents a vehicle front suspension response that has been measured using an accelerometer. It is seen that the measured rear suspension response is effectively just a time delayed version of the front suspension response and shifting the rear signal forward in time by the delay causes the signals to align rather well. The intention now is to simulate this rear suspension response without needing an accelerometer to measure it. Knowing the time delay between the two signals, I am trying to simulate the rear suspension response using a set of transfer functions that have been measured on the vehicle.

s1 is the front signal, s2 is the rear signal. s2sim is the rear signal that I want to simulate. h11, h12, h21 and h22 are the transfer functions used to try and simulate the rear signal. h11= a11/f1 where a11 is the acceleration measured at point 1 when a force is applied at point 1 and h12 = a12/f2 where a2 is acceleration measured at point 1 when force is applied at point 2 and so on. I am also looking at the phase information and making sure that the rear suspension response is replicated rather than the front suspension response.

S1=fft(s1);
Fop1= S1./(h11+(h12.*exp(-1j*w*timedelay)));
Fop2=Fop1.*exp(-1j*w*timedelay);
S2SIM=(h21.*Fop1)+(h22.*Fop2); % where Fop1 and Fop2 are the operational
                               % forces at the front and rear respectively
s2sim=ifft(S2SIM) % or real(ifft(S2SIM));

When I look at the phase information, it is nearly identical to the front suspension response rather than the rear and when i take the inverse fft and use the cross correlation to find out the time delay between s1 and s2sim in the time domain, it is seen that the time delay is 0 which means the signal hasn't been delayed so I am doing something wrong. Can anybody help me out here.


texts: Identifying the magnitude and impulse response from pole zero plot quickly

texts:
I have an exam next week and it's verty certain that a task of this kind will be there.

Are there some good tips how to match the right pole zero plot to the right responses? No proof is needed in the task.

At the moment I try to find two of the same kind and kinda guess which of the two same looking responses could be which.

What I know so far for the impulse response:

Amount of zeros = amount of poles -> dirac impulse at 0
Complex conjugated poles cause oscillation

For the magnitude response:

Zeros in the left half plane have low pass characteristics/highpass charactresistics in the right half plane
Poles outside the unit circle means unstable system

Are there quick ways without computing $\left\{\left|H(e^{j \Omega } )\right|\right\}$ ?

For A:
zero at z=1, and complex poles, so it has to be one of the oscillating impulse responses: 1 or 5. D has the poles closer together so it oscialltes less, so D is probably 1 and A is 5.

For B
Zeros are symmetric so magnitude response is probably symmetric.

For C
Amount of zeros = amount of poles, so dirac at 0 (2 or 3)

For D
See A

For E
No idea, should have a dirac at 0, but no oscillation?

Really unsure about the magnitude responses.


texts:
Hey Stack Overflow community,

I've been given a series of datasets that resembles a square pulse, and ideally looks like the following:

This particular plot is composed of many sets of data.

I'm wondering if there is a way to reduce the noise found in the data, and retrieve just the area under the square pulse with the assumption that "spikes" may appear before or on the actual pulse itself.

The method I've tried is to cut everything before and after the pulse, and look for spikes that may occur on the actual pulse by looking at how much each value deviates from the mode peak found. However, this seems less than adequate as my implementation requires going through each point multiple times generating quite a large time complexity, and the file sizes are pretty huge.

Thanks for helping out!


texts: I was reading a paper 'A Flexible Framework for Local Phase Coherence Computation' (paper URL) on using wavelets for local phase coherence. I just want to be sure. Is the gabor wavelet the same as the gabor function in the gabor transform (STFT)? If not, can we use a Windowed Fourier Transform instead?


texts:
I have the following transfer function
$$10000 \cdot\frac{(s-0.012)(s+1.05)(s-18)}{(s+0.22)^2(s+45)(s+1000)}$$
For which I am trying to tune a PID controller for. I'm using the pidTuner in MATLAB to do so. I wish to have the settle time be less than 3 seconds if possible, however, the best I have been able to do is 77 seconds with the settings $K_p = 0.049955$, $K_i = 0.0087525$ and $K_d=0$. At these values, going in any direction seems to increase the settling time more. Any help would be really great!


texts: What do different Cepstral or MFCC coefficients represent intuitively?

texts:
I understand the explanation for separating slow and rapid changing log spectral components but i need to understand:

Why lower coefficients have higher and mostly positive magnitudes?
Why Higher coefficients have  negative values more
often?

I understand the mathematical computation i.e. taking log of spectrum to make the components with different change rates, linearly seperable and then taking IFFT of logspectrum to separate these components.
However, I analyzed the MFCC stats from a speech dataset and can't see why slowly changing components (lower coefficients) have more magnitude while higher ones have more negative values and what do negative valued coefficients indicate.


texts:
I am trying to detect notes out of a musical sound input. My program can be run here. The problem is the patterns in the frequency response is repeating:

I suspect that echoing in the room is causing the higher frequency components. If I try a sound generator such as this, repetition is much less, but it is there. Try for example 165 Hz (note E3).
My primary question is that, how come sometimes repetitions have higher spikes in FFT than the actual signal (see the above image). This happens once in a while when human voice is tested. I know the actual signal is the 1st one because I test it in another software (Singer's Studio on iPhone).


texts: Why is this matrix invertible in the Kalman gain?

texts:
In the wikipedia article about Kalman filters, the well-known expression of the matrix of Kalman gains is given:
$$ \mathbf {K} _{k}=\mathbf {P} _{k\mid k-1}\mathbf {H} _{k}^{\text{T}}\mathbf {S} _{k}^{-1} $$
with
$$\mathbf{S}_k=\mathbf {H} _{k}\mathbf {P} _{k\mid k-1}\mathbf {H} _{k}^{\text{T}}+\mathbf {R} _{k}.$$

I understand that $\mathbf{R}_k$, as a covariance matrix, can be asked to be non-singular: it is reasonable to believe that no variance is zero. But this does not answer my question: why is $\mathbf{S}_k$ invertible?


texts: The phase of signal when it passes a frequency divider?

texts:
When we used the microwave frequency divider, such like: RF Bay FBS-N-14 or analog dividers in wiki.

The input signal is $\cos(\omega t+\phi)$, and the output signal is $\cos(\frac{1}{2}\omega t+\theta)$ when it passes the divide-by-2 prescaler.

what is the phase $\theta$ equal to? $\frac{\phi}{2}$ or $\frac{\phi}{2}+\pi$?


texts:
I am new to image compression topic and really appreciate it if someone could help me understand how compression rati works in this case. Assuming I have JPEG lossless coder that encodes the following image:

2 3 4 3
1 A B C

I have used predictive coding and Huffman to obtain the binary representations (codes) for pixels A, B and C:

Which are:

A=0 -> 0
B=3 -> 1101
C=0 -> 0

How can I calculate the compression ratio achieved for coding of the pixels A, B and C (jointly) assuming that the original image has pixel values ranging from 1 to 8?


texts: How to determine the central frequency from a vectorial spectrum analyzer?

texts:
I have taken measures of an antenna with a vectorial spectrum analyzer. I need to obtain the central frequency at which the antenna works. In order to do so, I have selected S11 (reflection) and I have plotted the "log magnitude" of the spectrum, I obtained something like this:

I know that the central frequency is obtained from one of the peaks shown on the plot. But in this case, I obtain three peaks. Which one should I take if they are pretty much the same? (no peak is considerably greater than the others)


texts:
I am trying to understand the recursive implementation of the DCT on the input signal (used for LMS filtering) according to the pictures below.

The pictures and formulas are taken from the paper: "Adaptive Inverse Control: A Signal Processing Approach", by Bernard Widrow and Eugene Walach

Block diagram of DCT structure:

Normal Matrix-Vector DCT-formula

Recursive DCT-formula:

What bothers me is how can this recursive DCT implementation yield the same result as the normal computation by using matrix-vector multiplication when the recursive formula uses the signal $x_{k-n}$?

But on the other hand the normal computation only uses signals up to $x_{k-n+1}$.


texts: How to find Integration of sinc function and then it's infinite summation?

texts:
I was doing this problem,i got answers of (a) and (b) by going following way,
first  i computed F.T. of x(t) as follows which in turn will be used to estimate it's bandwidth which will come as 2Hz and bandwidth of it's baseband version y(t) will come as 4Hz (if i am not wrong)
$x(t)=4sinc^2(2t-1) \xrightarrow{\mathcal{F}}2\Lambda \left(\dfrac{f}{2}\right)e^{-j\pi f}=X(f)$
such that,
$|X(f)|=2\left[1-\left|\frac{f}{2}\right|\right] \forall f \in[-2, 2]\\\\$
$\angle X(f)=-\pi f \ \ \forall f\in [-2, 2]$
and
$Y(f)=\Lambda \left(\dfrac{f-10}{2}\right)e^{-j\pi(f-10)}+\Lambda \left(\dfrac{f+10}{2}\right)e^{-j\pi(f+10)}$
But when Calculating Power of $z(t)$, which is I.F.T of frequency sampled version of baseband spectrum $X(f)$ ,in part $(c)$ i got stuck after proceeding as follows, and i want help regarding this (c) part only:
$P=\lim_{T \to \infty}\dfrac{1}{T}\displaystyle\int_{-T/2}^{T/2} |z(t)|^2 dt$
$=\lim_{T\to \infty}\dfrac{1}{T}\displaystyle\int_{-T/2}^{T/2}\left|8\displaystyle\sum_{k=-\infty}^{k=\infty}sinc^2(2t-1) \circledast \delta(t-2k)\right|^2 dt=64\displaystyle\sum_{k=-\infty}^{k=\infty}\lim_{T\to \infty}\dfrac{1}{T}\displaystyle\int_{-T/2}^{T/2}sinc^4 (2t-4k-1)dt$
(this is when i tried evaluating average power in time domain)
similarly,i also tried evaluating power in frequency domain by using following figure(bottomost right figure) but got no progress whatsoever

So, Please help me with the third part if possible.
Thanks for your response.


texts: Simulating Range Bearing Sensor with MATLAB with Gaussian Noise (Generating Gaussian Colored Random Vector)

texts:
I would like to simulate a sensor that provides range and direction of a beacon. This is for EKF localization, so the noise must be Gaussian (i.e. $\mathcal{N}(0, \sigma^{2})$. Also, I would like to compute the measurement noise matrix $Q_{k}$ which is
$$
Q_{k} = \begin{bmatrix} \sigma^{2}_{r} & 0 \\ 0 & \sigma^{2}_{\theta} \end{bmatrix}
$$

This is what I've done

sigma_range = 1;
sigma_angle = degtorad(5);
Q = [sigma_range^2,               0;
                 0, sigma_angle^2];

In the sensor, this is I've done

% add Gaussian noise with zero mean and some variance
Z(1) = Z(1) + sigma_range*randn(); % Z(1) the distance btw robot and beacon
Z(2) = Z(2) + sigma_angle*randn(); % Z(2) the angle btw robot and beacon

Is this correct?


texts: Alamouti code after or before IFFT block in OFDM

texts:
We are trying to implement Alamouti STBC on underlying ofdm signal. We were thinking as to where to write alamouti coder block.

I thought we could write it after the cyclic prefixer block as we get a symbol after the IFFT and cyclic prefixer. we could take two such symbols,send both of them at one time instant and then conjugate them and send them on the next time instant.

But then we saw a few block diagrams where alamouti coder is used  before IFFT block and after ofdm symbol mapper. when I searched why that's the case, there were reasons that said for OFDM to be effective, symbols should be assigned to subcarriers in frequency domain to mitigate multipath effects and not in time domain. The points that I don't understand are:

What does it mean by saying symbols in ofdm  are assigned in frequency domain on individual subcarriers? I can only see a symbol at the output of a ofdm mapper which is a complex number which comes from the constellation diagram depending on the modulation technique used. The output of the IFFT is also a complex number. So  how does an IFFT operation done on complex symbols mean that they are assigned in frequency. we could as well use a FFT instead of a IFFT.
What is the reason for alamouti code to be implemented before IFFT? Can it also be implemented the way I thought after cyclic prefixer block?


texts: How to prove: high-slope/more change is high freq part of signal

texts:
Generally, everyone knows that "the part of signal with high-slope is the high frequency part."

What is the definition or theory behind it?
How to prove this well-known thing?
perhaps, I miss something apparent.


texts: FIR filter (from zeros to filter coefficients)

texts:
Assume $F=1$ $(T=1)$. A 2nd order FIR filter with real coefficients has a zero in $H(e^{j\omega T})$ at $z_0=e^{j\omega /6}$. In addition, $|H(e^{j\pi T})| = 1$.

What are the coefficients of the FIR filter?

current progress

I'm having trouble understanding the $H(e^{j\omega T})$ function.
Tried to start from

$$y[n] = \sum_{i=0}^2 b_i z[n-2].$$

Does that fact that it has a zero in a point means I must apply zero-pole rules here?


texts:
I'm reading a paper on EMG analysis. The formulas are all clear to me, but the paper refers to the signal amplitude as "instantaneous" amplitude.

I know what instantaneous mean, but what does it mean in the context of signal processing?


texts:
I have a question about sammpling rate in MATLAB.
I read a .wav file using audioread(). It give me sampled data called my_data, and $88.2 \textrm{ kHz}$ as a samping rate, called Fs.

When I played this data with audioplayer(my_data, Fs), it works well and totally right sounds with original file.

What I tried to do is getting odd-indexed data like this :

my_data = [1 2 3 4 5 6 7 8 9 10]
my_odd_data = [1 3 5 7 9]

I thought that this procedure is like sampling my_data with half of Fs, but it is not true. I want to know why this procedure is not sampling with Fs/2.


texts: reconstruction filter - How does it actually work?

texts:
I'm trying to form my own understanding on the religious war around using 192kHz as a sampling rate for playback (the Internet seems to have a wealth of material on both sides). I'm struggling to understand how reconstruction filters work.

The Nyquist-Shannon sampling theorem, often cited by the anti-192 camp, basically states that a 44.1kHz sample rate is sufficient to reconstruct without loss a band-limited signal of 20kHz. However, when looking at the Whittaker–Shannon interpolation formula, it seems to me that an ideal reconstruction filter would need to have access to all samples, i.e., all past and all future samples.

I'm not an expert in analog audio, but I doubt that such a device can be built. I would assume that, at best, a delay can be introduced to "wait" for enough future samples to arrive, so that the contribution of unavailable future samples to the current output instant becomes negligible.

Can someone explain how practical reconstruction filters work and what are their tradeoffs? Is there a theoretically stricter limit to the Nyquist-Shannon theorem if only a window of samples is available or if reconstruction latency is unacceptable?


texts: Why do I have to scale filter coefficients during bandlimited sample rate conversion?

texts:
I am writing an example for how to perform sample rate conversion using windowed sinc and polyphase filters.

Here is the C++ code as a github gist.

So far, this works, and I think I have understood it correctly.

But, at line 184, you can see something curious. If the downsampling factor is higher than the upsampling factor (that is, we are overall downsampling the input by a fractional factor), I have to scale the filter coefficients, otherwise the output is too loud (its amplitude is increased). And I do not know why exactly that is. My current suspicion is that is necessary in order to account for the energy in the parts of the spectrum that I cut off if I overall downsample the signal (for example, when converting 48 kHz -> 40 kHz, this would be the energy of the upper 8 kHz that I cut off).

If I apply this scaling, it works, except for when I use very small filters (5 coefficients or less). With such small filters, the output signal gets dampened. I am unsure if this hints at a deeper problem, or if this is just to be expected with such small filters due to numerical errors.

Can anybody help out?

Oh, and there is also something unclear at line 169 regarding the parameter for sinc.


texts: Make a signal that fits another the best possible with a limitation in the 2nd derivative

texts:
Consider this step function:

The signal that "fits" this should look like the following (in green):

The corners are now smooth because the maximum second derivative allowed is not infinite anymore.

The idea is to find the most close to the original signal given the maximum second derivative allowed, minimizing the difference between the two signals (in orange):

If your system can predict things before it happens it could do something like this:

Minimizing the total difference even more. But also could be achieved if it works in non real time or has some kind of look ahead. But that will be not always possible or desirable.
This is not in the very core of my question, however is also a interesting option.

So, my question is: How can this be done? Is it possible (eg: exist more than one curve that minimizes error to the same value)?

Also to me looks like a filter, is there any filter which definition is something like this (hopefully a digital one)?


texts: Converting poles/zeros to differential/difference equation solutions

texts:
Does anyone have a reference handy on how to convert poles and zeros of a system to differential/difference equations.

Here is a quick draft of math, but I am not sure if it's at all correct. First, solving the characteristic equation gives the eigen values (equal to poles). Applying rudimentary knowledge of differential equations, the solution regarding only the poles should be:

$$\text {Poles Diffrential}: p(t)= \sum_{i=1}^{n_1} c_ie^{t\times \text{p}_i} $$

$$\text {Poles Difference}:p[n]= \sum_{i=1}^{n_1} c_i\text{p}_i^n $$

Where the $p_i$ is the value of ith pole and C is some constant, possibly scaled to the inverse value of the number of poles,  the index runs from 1 to $n_1$, the number of poles. I denote the solution that considers only poles/characteristic equation as $p(t)$. (To be exact, c is the input/initial value at time zero so we still have the denominator of characteristic equation as cz^0=c).

Next, zeros should be solved in respect to time/samples. This should be an easy task, considering that the input equation which zeros are solved from, directly outputs the signal. First solve the Y(z) (equivalent to the numerator of the transfer function):

$$Y[z]=\prod_{i=1}^{n_2}(z-q_i)$$

Where $q_i$ is the value of ith zero and $n_2$ the number of zeros. The function needs to still be converted to time doimain by converting each individual elements $az^{x}$ into $z[n(-x)]=a$. Things get more complicated when both elements are considered simultaneously. Unless my interpretation is incorrect, the input equation (used for zeros) acts as a signal that is fed into the recursive filter(poles portion). Thus the final interpretation can be had by simply convolving the equations (with $c_i=1$):

$$ \text {Zeros, Poles}: y[n] = p[n] * z[n] $$

Where the $p[n]$ and $z[n]$ from the preceding equations and $*$ is the convolution operator. To solve the differential equation variation, simply substitute $t$ in place of $n$ for zeros and the output.


texts:
I am getting 20-length mfcc of an audio signal. Each element(vector) is of size 5. I am using Librosa.
How to give custom size to the each element of mfcc (say 6)?


texts: Is there a closed form for the average energy of a QAM constellation?

texts:
I'm trying to come up with a formula for the average energy of the signals in a M-QAM constellation.
If $b$ is the number of bits represented by each point in the constellation, such that $M=2^b$. For $b$ even, I know the formula is:

$$ \bar{E} =  \frac{(M-1)d^2}{6} $$
where $d$ is the distance between two points in the constellation.

This can be easily verified using the average of all QAM constellation points generated by MATLAB, whose default value for $d$ is $2$:

b=4;
M=2^b;
sum(abs(qammod([0:M-1],M,0,'gray')).^2)/M

Is there a formula for the average energy in a QAM constellation with $b$ odd?


texts: Plotting Hilbert and Marginal Spectra in Python

texts:
I am using the PyEMD package for python 3.6. I want to create a 2d plot of hilbert spectrum (x: Time, y: Frequency, Colorbar: Amplitude) and the Marginal Spectrum (x: Frequency, y: Amplitude)

I can compute the instantaneous phase, frequency and amplitude of every IMFs:

import numpy as np
import pylab as plt
from scipy.signal import hilbert

from PyEMD import EMD

def instant_phase(imfs):
    """Extract analytical signal through Hilbert Transform."""
    analytic_signal = hilbert(imfs)  # Apply Hilbert transform to each row
    phase = np.unwrap(np.angle(analytic_signal))  # Compute angle between img and real
    return phase

def instant_amplitude(imfs):
    """Extract analytical signal through Hilbert Transform."""
    analytic_signal = hilbert(imfs)  # Apply Hilbert transform to each row
    analytic_amplitude = np.abs(analytic_signal) ## Compute amplitude
    return analytic_amplitude

# Define signal
t = np.linspace(0, 1, 200)
dt = t[1]-t[0]

sin = lambda x, p: np.sin(2*np.pi*x*t+p)
S = 3*sin(18, 0.2)*(t-0.2)**2
S += 5*sin(11, 2.7)
S += 3*sin(14, 1.6)
S += 1*np.sin(4*2*np.pi*(t-0.8)**2)
S += t**2.1 -t

# Compute IMFs with EMD
emd = EMD()
imfs = emd(S, t)

# Extract instantaneous phases, frequencies and amplitude using Hilbert transform
instant_phases = instant_phase(imfs)
instant_freqs = np.diff(instant_phases)/(2*np.pi*dt)
instant_amplitude = instant_amplitude(imfs)

Now I am stuck and don't know how to compute and plot these formulas:

Wikipedia: Hilbert Spectrum

Can someone help me out?


texts:
I'm after some advice on a problem that I am facing regarding signal processing.

I am processing a voice clip (.wav) and this voice clip contains noise, and I need to remove the noise and just concentrate on the Phones so I can analise them.

What I have thought about, and, what I have researched is, if I use zero-crossing as well as the total energy I can remove the blocks that do not contain sufficient energy or frequency.

The problem at the minute is that I have a 1D vector, acting as a 2D vector and let's assume the following data is used:

v1 = {128, 128, 128, 127}

Now my current algorithm squares each of these numbers and then sums the entire vector up. This, would therefore not be able to give me an accurate identification of noisy signals. So should I therefore split the signal into blocks (In that I get a 2D vector)? e.g.

B1 = {128, 128};
B2 = {128, 127};

e1 = {sum( (128^2) (128^2) )};
e2 = {sum( (128*2) (127^2) )};

Therefore, this would mean I am given 2 (In this example) energies which I can then check to see if they match the threshold.

Last question, do you think the methods of zero-crossing AND energy are good an accurate way of determining noise and removing it?

Hope someone can help me, thanks :)!


texts: Amplitude Estimating Using a Windowed DFT

texts:
Let's say we want to estimate the amplitude A of a mono-frequent signal using a windowed DFT. The frequency of the signal is unknown, and the frequency resolution of the DFT is limited, thus it cannot be expected that there is a DFT bin that matches this frequency directly. Therefore, A is estimated from the amplitude of the next closest DFT bin. Now the question is,

should the frequency response of the window exhibit a narrow or a wide mainlobe to obtain a good estimate, and why?


texts: Counting human movement of exercise using ROI or pre-set passing line with OpenCV

texts:
I'm currently doing my final year Degree in computing studies.

The project I'm working on is a video of a webcam or computer webcam to have 2 pre-defined ROI (Region of Interest) on top and bottom of the video. It is to detect a person when he does a simple push up whether he reaches the area of both top and bottom ROI by detecting the person's head when he enters the ROI. I totally have very less experience in programming and really need guide and help.

During my research I've come across a journal/research at: Pattern Recognition Letters 33 (2012) 438–445. I've downloaded and gone thru it but can't seem to find its code or truly understand it. I got access with my University account at sciencedirect.com.

2nd idea I had was after viewing a video in Youtube: Traffic Counting based on OpenCV... I was very interested by that video. I was thinking if I could set 2 lines like he did and detect simple human motion up and down doing a simple accurate pushup exercise I would be so much happy but I can't seem to get contact or ask for any information.

I really need advice and help because I am really lost.


texts: Failed to implement Goertzel algorithm in Python

texts:
After some questioning on stackoverflow, I tried to implement a Goertzel algorithm in Python. But it doesn't work : https://gist.github.com/4128537

import math

def goertzel(samples, sample_rate, f_start, f_end):
    """
    Implementation of the Goertzel algorithm, useful for calculating individual
    terms of a discrete Fourier transform.
    """
    window_size = len(samples)
    f_step = sample_rate / float(window_size)
    # Calculate which DFT bins we'll have to compute
    k_start = int(math.ceil(f_start / f_step))
    k_end = int(math.floor(f_end / f_step))

    if k_end > window_size - 1: raise ValueError('frequency out of range %s' % k_end)

    # For all the bins between `f_start` and `f_end`, calculate the DFT
    # term
    n_range = range(0, window_size)
    freqs = []
    results = []
    for k in range(k_start, k_end + 1):
        # Bin frequency and coefficients for the computation
        f = k * f_step
        w_real = 2.0 * math.cos(2.0 * math.pi * f)
        w_imag = math.sin(2.0 * math.pi * f)

        # Doing the calculation on the whole sample
        d1, d2 = 0.0, 0.0
        for n in n_range:
            y  = samples[n] + w_real * d1 - d2
            d2, d1 = d1, y

        # Storing results `(real part, imag part, power)`
        results.append((
            0.5 * w_real * d1 - d2, w_imag * d1,
            d2**2 + d1**2 - 2 * w_real * d1 * d2)
        )
        freqs.append(f)
    return freqs, results

if __name__ == '__main__':
    # quick test
    import numpy as np
    import pylab

    t = np.linspace(0, 1, 44100)
    sine_wave = np.sin(2*np.pi*441*t)[:1024]
    freqs, results = goertzel(sine_wave, 44100, 0, 22049)
    print np.array(results)
    pylab.plot(freqs, np.array(results)[:,2])
    pylab.show()

I am a beginner in this topic, so I have no idea what could be wrong in there. Any advice would be welcome.

EDIT

Here is what I get when plotting the power ... as you can notice, the frequency 440 which should appear is not there :


texts: How to convert magnitude to dB in matlab

texts:
I exported a diagram from FDTD simulation, and now I want to change the magnitude to dB, here is the exported codes:
ax1 = axes;

plot(lum.x0,lum.y0)

set(ax1, 'XLim', [1.49806e-06 1.60194e-06])

set(ax1, 'YLim', [-0.0232806 1.22328])

set(ax1,'XGrid', 'on')

set(ax1,'YGrid', 'on')

and the ploted spectrum is shown below

I just wondering how to convert the magnitude to dB?
I know there is a function called mag2dB, but i just dont know how to use it.
any help would be very appreciate.


texts: What are the cosine functions in JPEG's DCTII table?

texts:
As I understand it from this video, in a JPEG image an 8x8-pixel block is made up of weighted cosine waves, calculated using DCTII. There are lots of visualizations of these waves, such as this one from the Wikipedia article on JPEG compression:

I'd like to be able to visualize each of the waves for a given 8x8-pixel block in an image, weighted with values calculated by DCTII.

Where can I find (or how can I calculate) the 64 cosine functions used in JPEG compression?


texts: I have a long audio file (in wav) and a short audio file (in wav) that has been copied from the long file. Now I need to find the correct location where the short signal has been copied from. I tried to compare the sample data but for some reason (possibly resampling or introduction of short fades) it doesn't match up. What is a good and fast way to solve this?


texts: Finding radius r of the overlappable sphere(s) in 3D image

texts:
My current problem:

I have an input 3D binary image (a 3D matrix that has only 0 and 1) that consists of random numbers of sphere with radius r.
We do not know how many spheres are there in the image.
All spheres have the same radius r, but we do not know radius r.
Spheres are all over the place in the image and can overlap each other.
example image is given below.

My requirement:

what is the radius r?

Currently, I simply just flatten the image to get rid of the z axis and perform edge detection and I am trying Hough Transform using: http://rsbweb.nih.gov/ij/plugins/hough-circles.html

However, with Hough Transform, I see that the variables minimum radius, maximum radius, and number of circles have to be specified. I have tried a few attempts below:

Given the correct parameters, Hough Transform can detect the circles just fine. But in the real application, I do not know how many spheres are there, and making the program attempting to guess minimum and maximum radius seems not feasible. Are there other ways to accomplish this?

Cross-link: https://math.stackexchange.com/questions/118815/finding-radius-r-of-the-overlappable-spheres-in-3d-image


texts: Single Pixel Camera - Compressive Sensing

texts:
I work with a setup that measures high frequency fluctuations in light using a photodiode. We  steer the light over a sample as we measure these fluctuations.

I am familiar with compressive sensing where we only sample a few points to reconstruct the full signal. The sample is of course done in the same dimensions as the signal - for audio, the sampling is across time for instance. For an image, we randomly subsample in both x and y coordinates.

I came across a single pixel camera which employs a Digital Micromirror Device to project a pattern onto a scene which is then imaged with, again, a photodiode. Most of the descriptions I've read say that the each pattern would in effect sum the light from the scene and repeated scanning with random patterns can yield images with similar quality as raster scanning.

How are these images reconstructed? If we sum the light, we lose where the individual sums of that light came from. How does repeated acqusition solve this?


texts: Why is In-phase Quadrature sampling not used to record and store digital audio?

texts:
I/Q is a standard method to sample radio frequency and its use is even seen in some non radio signals like sonar. But it seems that it is not used to sample acoustic frequencies.

What property of sound makes plain in-phase PCM or PDM adequate for recording and storing audio–even for high-fidelity recording?

I have an intuition as to the answer, but I think satisfying explanation would be useful to those researching this area of signal processing.


texts: Maximum likelihood estimator of active time delay and passive time delay

texts:
A typical time delay estimation problem has the model:
$$
\begin{align}
     x_1(t) &= s(t)+ n_1(t) \\
    x_2(t) &= a s(t-D) + n_2(t)
\end{align}
$$
Where $n_1$ and $n_2$ are considered to be Gaussian noise.
An active time delay estimation requires $n_1$ to be $0$ and $a=1$.

The maximum likelihood estimator, in this case, is the matched filter or the basic cross-correlation method for finding delay. In passive case, $n_1$ is not equal to 0, and the maximum likelihood estimator is HT processor in Generalized Cross Correlation (GCC).

Can someone provide me with an intuition why in active time delay estimation case, the matched filter performs favorably better than GCC, even if I know the spectra of my signal (i.e. the spectral estimates are accurate enough)?

I would expect both GCC and simple CC to give almost equal results. And also, my results vary when estimating time delay by cross correlation in time and frequency domain. Is there any explanation behind this?


texts:
Are there any established techniques for estimating the envelope of a swept filter? I tried a first order piece-wise polynomial approximation using the SHGO global minimizer algorithm (Python code below), but it fails. Setting the points by hand approximately succeeds.
The purpose for the estimation is to obtain a synthesizer's filter decay envelope, after the filter is known. The input is thereby one of the typical waveforms a synthesizer can generate, in the code a sawtooth waveform was chosen.
import numpy
import scipy.optimize
import scipy.signal

#2nd order low pass filter
def filt(X,F,X0=0,X1=0,Y0=0,Y1=0,resonance=1):
    X=numpy.append(X,[X1,X0])
    Y=numpy.zeros(len(F))
    Y=numpy.append(Y,[Y1,Y0])
    for i in range(len(F)):
        wd=F[i]*2*numpy.pi
        T=1/88200
        d=resonance
        C=numpy.tan(wd*T/2)

        D=(1+d*C+C**2)
        a1=2*(C**2-1)/D
        a2=(1-d*C+C**2)/D

        b0=C**2/D
        b1=2*b0
        b2=b0

        Y[i]=-a1*Y[i-1]-a2*Y[i-2]+b0*X[i]+b1*X[i-1]+b2*X[i-2]
    return Y[0:len(Y)-2]

#X = filter input, Y = Output.
freq=600
t=numpy.linspace(0, 1, 88200)
X=scipy.signal.sawtooth(2 * numpy.pi * freq * t)
F=freq+5*freq*numpy.e**(-t*40)      #Envelope
Y=filt(X,F)
Y+=numpy.random.normal(0,0.075,len(Y))

#Estimation (Result in Fest):
def PolMin(X,P,I,X0,X1,Y0,Y1):
    F=numpy.zeros(I[1]-I[0])
    F=numpy.polyval(P,range(I[1]-I[0]))
    Y=filt(X[I[0]:I[1]],F,X0,X1,Y0,Y1)
    return Y

def Target(P,I,X0,X1,Y0,Y1):
    return numpy.sum((PolMin(X,P,I,X0,X1,Y0,Y1)-Y[I[0]:I[1]])**2)

chunkl=200
length=25000
chunks=int(25000/200)
Results = numpy.zeros([chunks,2])
B=[(-2,0),(0,10000)]
Y0=0
Y1=0
X0=0
X1=0
Z=numpy.zeros(length+1)
Fest=numpy.zeros(length+1)
for i in range(chunks):
    chunkst=i*chunkl
    chunke=chunkl*(i+1)+1
    f=lambda x: Target(x,(chunkst,chunke),X0,X1,Y0,Y1)
    Results[i,:]=scipy.optimize.shgo(f,B,iters=4).x
    endpoint=numpy.polyval(Results[i,:],chunkl*(i+1))
    B=[(-1,1),(endpoint,endpoint)]
    Fest[chunkst:chunke]=numpy.polyval(Results[i,:],range(chunkl+1))
    Z[chunkst:chunke]=filt(X[chunkst:chunke],Fest[chunkst:chunke],X0,X1,Y0,Y1)
    Y0=Z[chunke-2]
    Y1=Z[chunke-3]
    X0=X[chunke-2]
    X1=X[chunke-3]


texts: Improving credit-card-like image scan for OCR

texts:
OCR works best for black on white or binary images, where letters have sharp contours, no noise, and solid filling, i.e. they are nothing like this:

For the curious it's a part of a document number that is embossed.
Applying unsharp mask, thresholding, closing and median operations already improves it a bit, but it still cannot be processed properly by OCR.

What other operations could bring the image to the preferable state like this manually edited one:


texts:
I have the dynamics of a 2nd order system, mass-spring-damper for example, in the transfer-function format.
For the analysis that I am doing, I am calculating the frequency response of the system by replacing the $s$ with $j\omega$ for a given range of $\omega$ of interest.
Now, I need to examine in the behavior of the system when noise is added. How does one go about adding noise to the frequency response? Is it valid to generate some random noise (constant mean and variance) divided into a real and imaginary part and add it to the frequency response?


texts:
Hello can someone explain the semantical meaning of the structure tensor and the Hessian matrix. I am aware of how it is calculated, but i find it difficult to comprehend what they describe in the image and how they can be used to detect corners (Harris) and detect Blobs (SURF).

Any help is greatly appreciated, thank you in advance.


texts: Wiener filtering/deconvolution for non-stationary noise

texts:
Consider a stationary discrete-time random process $x[k]$ which undergoes low pass filtering by a filter with impulse response $h[k]$ and is subject to additive, temporally uncorrelated noise $n[k]$ of time-varying power $\sigma^2_n[k]$:
$$
y[k] = \sum^\infty_{m=0} h[m] x[k-m] + n[k].\qquad
$$
The time-variation in noise power is described by:
$$
 E[n^2[k]]=\sigma^2_n[k]= Ae^{k/B} \qquad B>>1
$$
such that the noise power increases (slowly) with time.
Specification of length $L$ FIR filter coefficients of a Wiener deconvolver seeking to estimate $x[k]$ from $y[k]$ is straight-forward for stationary noise.
Things are more complicated for non-stationary noise.
Strictly speaking non-stationary noise requires a different filter for every sample instant $k$ which is prohibitively complex.  In my case, the noise changes slowly enough for it to be modeled as locally stationary over significantly more than $L$ samples.  I could always just update the filter coefficients every $NL$ samples to reflect a new noise power assumption, but I'm concerned about abrupt transient responses due to the sudden jump in deconvolver coefficients.  Something more frequent but more gradual sounds better.
Is anyone aware of a smart recursive form of the Wiener filter for cases like this (especially where the noise power itself follows a recursion: $\sigma^2_n[k]=e^{1/B}\sigma^2_n[k-1]$)?


texts:
I have come up with a new descriptor(floating point). I want to evaluate this descriptor against SIFT, GLOH, etc etc. I am using the K. Mikolajczyk, C. Schmid, approach.  The method is described in the link http://www.robots.ox.ac.uk/~vgg/research/affine/descriptors.html. They have provided the binaries and the data set. I extract the patches(41x41) and apply my descriptor on these patches. But the problem is, these patches are not rotational normalized.

SIFT achieves rotation in-variance by rotating the patch in the direction of the dominant angle of the patch. My query is ,how to rotate the patch or achieve rotation normalization? I want the original patch(41x41) and the rotated patch to have the same dimension (41x41). For example, If i use the matlab function imrotate and rotate the original patch by 60 degrees. The rotated patch has a dimension bigger than (41x41). And there are blank spaces in the border of the rotated image. I think applying my descriptor on these rotated and enlarged patches gives wrong results.

M question is, How do i rotate the patches or rotational normalize the patch(Such that the original patch and the rotated patch have the same size 41x41) to apply my descriptor?


texts: Converting transfer function that is a sum of unusual rational polynomials to finite difference equation

texts:
I have the following rather exotic transfer function:

$$
H(z) = cz^{-m} + \frac{b_0 z^{-1} + b_1z^{-2} + \dots + b_{2m}z^{-2m}}{1 + a z^{-1}} + \frac{q_0 z^{-1} + q_1z^{-2} + \dots + q_{2m}z^{-2m}}{1 + p z^{-1}}
$$

Where $c$, $b_0,\ldots,b_{2m}$,  $a$, $q_0,\ldots, q_{2m}$, $p$ are real valued constants and $m > 1$.

How do I go about converting this to a finite difference equation via the inverse $\mathcal Z$-transform?

I thought I had got to grips with the inverse $\mathcal Z$-transform, but I can't seem to convert $H(z)$ to a finite difference equation that produces the right output.

My attempt:

$$\begin{align}
y[n] &= c x[n - m] &+& (b_0 x[n - 1] + \dots + b_{2m} x[n - 2m]) \\
&&-& a y[n - 1]\\
&&+& (q_0 x[n - 1] + \dots + q_{2m} x[n - 2m])\\
&&-& p y[n - 1]
\end{align}$$

This doesn't seem to be correct. When I test this its magnitude response is not even close to the original transfer function.

My assumption is that the three terms $cz^{-m}$, $\frac{b_0 z^{-1} + b_1z^{-2} + \ldots + b_{2m}z^{-2m}}{1 + a z^{-1}}$, $\frac{q_0 z^{-1} + q_1z^{-2} + \ldots + q_{2m}z^{-2m}}{1 + p z^{-1}}$ can be transformed separately.


texts: How the weight $ω_{i,t}$ is initialized in Background subtraction GMM (Gaussian Mixture Model)?

texts: I read Chris Stauffer paper "Adaptive background mixture models for real-time tracking" about Gaussian Mixture Model for background subtraction, but they're not explained very detail about $ω_{i,t}$ value, what I want to ask is how the value initialized and what the range of it?


texts: Is it possible to do further beamforming on already beamformed data?

texts:
I've a circular array of eight FM antennas. I'm doing digital beamforming to form eight beams with the following beam patterns:

The objective is to reduce my side lobes even further, and to get narrower main beam. What I'm interested in to know if it is possible to take two adjacent beams, and do another step of beamforming, i.e., considering them as if coming from two antenna array with the beam pattern shown in the figure. I don't know how to incorporate the beam pattern into the standard set of equations for beamforming, as I think an omni-directional antenna is assumed. Any directions, suggestions are welcome.


texts: Is it possible to get the velocity vector for tangential movement using radar?

texts:
When an object appears towards the radar sensor or goes away from the radar sensor, there is a change in the phase shift of the frequency and from that, the velocity can be inferred.

So, I have two questions.

Is it also possible to find out the velocity "vector" ( ie which direction the object is moving to and not just the magnitude of the velocity of the object )
If the object is neither moving towards the sensor nor away from the sensor ( moving tagentially ) , will there be still a phase shift on the radiated pulse?


texts: Estimating the time offset (delay) between two audio signals in real-time

texts:
I'm looking for an algorithm which can accurately estimate the time-offset/delay between two audio signals in real-time. One audio signal is a loudspeaker signal and the other signal is the microphone signal.

Audio is processed in blocks of 10ms (sample frequency 16kHz).

There is light office noise embedded in the signals
The algorithm must be able to estimate delays up to 500ms.
The delay is usually constant for a period of time (with some variability -10ms to 10ms), but occasionally the delay "jumps". Those sudden jumps can be in the order of 50 to 80 ms. The algorithm must be able to handle that.

I've read that GCC PHAT is one algorithm which can be used, but I'm not sure how it would work in a real-time environment.


texts: Noise Removal from an Image Using OpenCV (Comparison to Neat Image)

texts:
I tried removing noise from the image shown below using Median Blur in OpenCV. But i'm not able to remove the colour noise completely as it is done in Neat Image. Any suggestions.?
 1. Original Input Image

 Median Blur Output
 Neat Image Output


texts:
I found this equality in the middle of my DFT calculation. Can anyone help me to proof this.

$ \frac{1}{N} \sum\limits_{k=0}^{N-1} \left(-W_N^n \right)^k \ = \ \delta[n-0.5 N] $

where

$ W_N = e^{j 2 \pi/N} $


texts: What's the point of MATLAB's scaling array for secord-order-section filters?

texts:
I'm currently working on a project where I'm implementing different filter types (e.g. iirs - direct form 1, 2, transposed 1, 2, etc.). I just recently started attempting to implement a second order section filter. This filter is made up of a cascade of biquad filters each implemented as transposed direct form 2 stages.

After implementing this, I wanted to verify that my results matched those in MATLAB. I opened up the filter design toolkit and generated an 8th order SOS filter. When I exported this filter, it exports both a SOS matrix along with a scaling array G.

So my question: Why does MATLAB provide this scaling array G? Shouldn't the SOS matrix be fully capable of describing the filter by itself? Why does MATLAB decide to factor out these gains?

Is this just a courtesy of MATLAB, such that the user can 'fine tune' each filter stage gain to determine the best option?


texts: What is localizability in Computer Vision?

texts:
Please consider the following excerpt (from Rodrigo R, Zouqi M, Chen Z, Samarabandu J.: Robust and efficient feature tracking for indoor navigation):

  Robust feature tracking is a requirement for many computer vision tasks such as indoor robot navigation. However, indoor scenes are characterized by poorly localizable features.

I am trying to understand what localizability is, and would appreciate help.  (PS: I know there is a dsp tag called "local-features", but the tag currently has no definition.)


texts: Using NMF for guitar pitch detection using fixed pitch template

texts:
I've seen this guide showing a basic NMF of a spectrogram.

My aim is detect guitar pitch from the decomposition. I want to do this by explicitly create a fixed template of pitches and construct the W matrix from these.

I've seen some papers giving an overview of this, but I am not sure how it's practically done.

Do I record every note on the guitar, then do the FFT for each, take its spectrum, and paste it inside a matrix?

I guess my question is basically: how do I create the fixed W matrix, given a set of recordings with each pitch.

This paper says:

  The whole note sample k is first processed in a short-time sound
  representation supposed to be non-negative and approximatively
  additive (e.g. a shorttime magnitude spectrum). The representations
  are stacked in a matrix $V^{(k)}$ where each column $v^{(k)}_{j}$ is the sound
  representation of the j-th time-frame.

So, does that mean we compute the STFT of the recorded sample? Then we use the resulting spectrogram as $V^{(k)}$?

Also, does the duration of the sample matter? If there is absolute silence in the beginning of the sample, will that affect $v^{k}_{j}$?


texts:
I have been doing some research into modelling acoustic spaces, and (with help from this site) I've managed to take the frequency response of fibrous layer model, convert it to a reflection response as a target for an IIR, and then find the polynomial coefficients for it using FDLS.  All good so far.

I performed the usual BIBO stability test of finding the magnitude of the IIR denominator roots, and with a 4th order filter everything was fine (<=1.0).  But increasing the order, even by just one, then it became unstable.

My question is: What does the order actually represent response-wise?  I always assumed it determined the sharpest change in direction of the frequency response i.e. the 'jagged-ness' or complexity of the response curve.  If that's true, then my filter is unstable, but the few poles I've calculated just happen to be inside the unit circle.

But is that true?


texts: Linear Convolution in DSP and Hann window

texts:
My understanding is that when dealing with discrete finite sequences, convolution done thru FFT is circular. To obtain linear convolution, one would pad the input with zeros up some appropriate length. Then take FFT of two padded sequences, element-by-element multiplication, then IFFT.

For very long signal, to do Linear convolution, one might use the Overlap-Add (OLA) method, which breaks the long signal into segments of equal length, then do the linear convolution as above and sum all the overlapping section together to get the final output.

My question is that isn't that when one breaks a long signal into parts and to use FFT on it later, one should apply a Hann window to minimize Gibbs phenomena due to discontinuities that might occur at two end points ?? So far all the documents about OLA that I come across has not shown use of Hann window in conjuction with OLA. I wonder if it's somehow not needed or all the OLA document just want to illustrate one idea clearly, so as not to confuse reader with addition of Hann window. I understand that if we use both OLA and Hann, it will be more complicated. Please enlighten. Thank you.


texts: Does Wav format supports writing complex values

texts:
I have a sine tone, doing some math, and finally gets an harmonically accumulated sine tone in complex numbers. I do want to preserve the phase information of imaginary part on the other end, which will continue processing the complex signal. So,

Is there any way i can safely pack the complex number into any of standard file formats, such as wav ?

Note : Currently I am  writing this signal's byte values to a .bin file


texts: How Is the Gaussian Kernel Related to the Euclidean Distance of the Neighborhoods in the Non Local Means Denoising Algorithm?

texts:
I am reading the paper A non-local algorithm for image denoising that describes the original non-local means algorithm. They define (p. 3 of the pdf) the distance between two square fixed neighbourhoods of two pixels as

\begin{align}
\|v(\mathcal{N}_i) - v(\mathcal{N}_j) \|_{2, a}^2
\end{align}

where $a> 0$ is the standard deviation of the Guassian kernel. However, I don't understand this notation. How is the Gaussian kernel used in the formulation? How is it related to the Euclidean distance?

I would appreciate if someone more familiar with image processing and computer vision and all this notation can give me an explanation.


texts:
Is there any better measurement of error between two images than just taking pixel value differences?

The problem I am tackling is stitching images containing "sensitive" objects, such as moving people. Here are two images I would like to blend:

I compute error map using inter-pixel differences:

The image is then considered a graph where pixels are nodes and edges are connecting neighboring pixels. Edge weights are differences between pixels taken from the error map.

Any Max-Flow/Min-Cut algorithm will find an optimal seam which goes through minimum error path (minimum cut - highlighted in red):

As you can see, the seam have successfully avoided girl's hand, but still goes through the guy's face.

The problem is in the error map which assigns too low error for his face as he just moved his head and the pixel differences of his skin are smaller than ripples on water and other less important features.

I heard about saliency map, but not sure if it is the right approach.


texts:
I need to study a wave file for every 100 Hz. However, I dont know how can I design filer for each 100 Hz. Let's say 0-100 Hz, 100-200 Hz, ... , ((fs/2)-100) _  (fs/2) Hz.

Is it any MATLAB code for it ?
How many samples are needed for bandwidth of 100  Hz ? I maen how should calculate it ?
Please guide me.
Thanks.


texts: NLMS Echo Cancellation: How do we estimate when the time-step at which the Far-End Echo will be generated?

texts:
I'm new to this domain. From my understanding, we do the following in our algo:

Error = Desired - Expected Output
That is,
E = (Near End voice + Far End Echo) - (Far End voice * Adaptive Weights)

Here, I assume the Far End voice is a block from a previous timestamp and the Near End voice and Far End Echo are blocks from the latest timestamp.

So for the application of this formula, how do we estimate after how many timestamps the Echo will be generated? That is, how do I know which timestep block should I use as the Far End voice to compare with the latest timestep input to update my Adaptive Filter weights?

Thanks in advance.


texts: Determining The Sparsity or Density of a Chroma Vector

texts:
I have a range of songs where I have took an 'unfolded' chromagram computed from a series of short-term windows.  Specifically, for one song, I have one vector which contains the amplitudes of the chroma's, and a second vector which contains the corresponding MIDI notes.  For now, I can plot a histogram of the amplitudes w.r.t the MIDI notes that gives me a visual idea of the range of notes used.

Here are two songs graphed:

As we can see, the 'classical' example is in my opinion more sparse than its neighbouring graph.  I've been trying to figure out a way as to how to measure for something like this, since I need to fold such a feature vector into 1 useful statistic.


texts: How can I plot a 3D graph of a given Laplace Transform of a function?

texts:
Let's say I have a function called $f(t)$ in time domain as:

$$f(t) = \exp(-3t)\cos(5t)$$

And the Laplace transform of this function call it $F(s)$ becomes:

$$F(s)=\frac{(s + 3)}{(s + 3)^2 + 25}$$

I want to plot the 3D plot of $\lvert F(s)\rvert$ as a surface above the $s$-plane.

I couldn't find any script. I have the version R2014a of MATLAB.
How is this done in MATLAB? Something similar to this plot:


texts: Understanding SURF point-of-interest localization

texts:
I've recently read Herbert Bay's paper of SURF. From which I understand:

SURF calculates a table from each image pixel's Hessian matrix's
determinant, instead of SIFT's Difference-of-Gaussian table. Then
SURF finds local maximal (not extremal) by comparing each determinant
to its 8 neighbors in the same scale, and 18 neighbors from 2 nearby
scales. Finally, SURF interpolates the maximal point to find the
exact position of the point of interest.
In the detecting part, SURF uses 9×9 approximations of second-order
Gaussian derivatives, combine with an integral image to fast
calculate the determinant of each image point's Hessian matrix. The
calculation for each derivative of the Hessian matrix is done by
placing the correct 9×9 filter center point to an image pixel x,
and simply multiplying corresponding points of the filter and the 9×9
sub-image whose point x is the center (points that are off image
equal to 0), then calculating the sum of all those multiplications.

Did I take something wrong? Any help is appreciated.


texts:
Consider a QAM system designed to transmit over a bandwidth of $3 \ \mathrm{kHz}$. The channel's power constraint imposes a maximum $\mathrm{SNR}$ of $30 \ \mathrm{dB}$. The system can tolerate a probability of error of $10^{−6}$. I want to determine the maximum throughput of the system in bits per second.

Let's call the throughput $T$. If $M$ is the constellation size and $"$ is the bandwidth, then I know that $T=MW$. I've calculated the value of $M$ using theformula $$M = \log_2 \left( 1 - \frac{1.5\mathrm{SNR}}{\ln(p_e)}\right)$$ Then I approximated $M$ to the nearest power of $2$, which is $8$.
Therefore, $T = 24000 \ \mathrm{\frac{bits}{sec}}$
but it is a wrong answer. Can anyone tell me what I have done wrong?


texts: Comparing similarity of two simultaneously recorded EEGs in freq domain

texts:
I have two EEG signals that are simultaneously recorded during the night (sleep). Imagine one device recorded from channel F4 and the other from Fp2. Now I would like to check the similarity of these signals, so, the first step is to sync them which has been done successfully using cross-correlation :)
I can also report some similarity metrics in time-domain, e.g. Pearson correlation per window (e.g. each 30 s) to quantify the similarity of signals.
The problem arises when I want to compare (or let's say quantify) them in the freq domain. When I plot e.g. spectrogram or PSD we can see their similarity visually; however, I don't know how to quantify this similarity in the frequency domain. By vision the spectrograms (or PSDs) are quite similar but I guess it is not a powerful reason without reporting any similarity metric as a number. Is there any metric by which I can report similarity of PSDs or spectrograms window-by-window and based on each frequency bin (Delta, Theta, etc)?


texts: Noise removal using adaptive noise cancellation algorithms in real time systems

texts:
I am newbie in signal processing. I want to remove noise from real time audio signal using adaptive noise cancellation(LMS&RLS) algorithms. I don't know how to choose reference signal(noise signal) for the algorithm from the incoming noise corrupted signal.  So can anyone help me to select the noise signal for the reference of ANC algorithms?
Thanks in advance


texts: Trouble with fixed-point arithmetic 2nd order IIR implementation in c

texts:
DSP newcomer here!

I am tinkering around with a TI DSP and am trying to implement a second order IIR filter in C. Input is 16 bit 2's complement, as is the output, the accumulator is 32 bit wide.
I have tried several structures: Direct form I, Transposed direct form II, Gold & Rader structure. However, I can't seem to get anything to work and I am not sure whether it is a typo in the code of all of my filter implementations or whether it is because of me misunderstanding the principle of fixed-point arithmetic. I assume - or better: hope - that it's the same mistake that is affecting all implementations.

Here's my code for the transposed df II implementation:

// state variables
int w[3] = {0};

short filter_dfii(short value) {
    w[0] = ((value * sos_filter_coeffs_b[0] + w[1]) >> 15 ) & 0xffff;
    w[1] = value * sos_filter_coeffs_b[1]
       - w[0] * sos_filter_coeffs_a[1]
       + w[2];
    w[2] = value * sos_filter_coeffs_b[2]
       - w[0] * sos_filter_coeffs_a[2];

    return (short) w[0];
}

My df I implementation looks like this:

int last_value[2] = {0};

short filter_dfi(short value) {
    w[0] =
        ((value * sos_filter_coeffs_b[0]
        + last_value[0] * sos_filter_coeffs_b[1]
        + last_value[1] * sos_filter_coeffs_b[2]
        - w[1] * sos_filter_coeffs_a[1]
        - w[2] * sos_filter_coeffs_a[2])
        >> 15) & 0xffff;
    last_value[1] = last_value[0];
    last_value[0] = value;
    w[2] = w[1];
    w[1] = w[0];
    return (short) w[0];
}

The filter coefficients are:

int sos_filter_coefficients_a = {32767, -51150, 21015};
int sos_filter_coefficients_b = {658, 1316, 658};

How did I get those coefficients? I utilized the Matlab 'butter' function and specified a 2nd order Butterworth lowpass with a cutoff frequency of 0.05*fs or 2400 kHz in my specific case. Then I multiplied each coefficient with 32767 and rounded it to an integer value. If I call Matlab freqz function with those quantized coefficients, I get a bode plot with the expected result, so I believe quantization error is not the issue here.

However, the result is not a lowpass but just garbage. Unfortunately, I lack the proper testing equipment to describe the result scientifically, but it sounds like digital oscillation that is slightly modulated by the actual input signal and is present even when there's no input at all.

So, in my opinion, this can't be an unstable filter since the poles are well within the unit circle. It can't be overflow limit cycle either because then it shouldn't be happening when there's only low input level or no input at all.. right?

If I alter the coefficients to

int sos_filter_coefficients_a = {32767, 0, 0};
int sos_filter_coefficients_b = {32767, 0, 0};

input gets returned basically unprocessed as expected, so my issues most likely aren't caused by any external factors - there must be something wrong with my coefficients, with the function itself, or both.

Any help is much appreciated!


texts: Do Nyquist samples have the same sum as the integral of the signal?

texts:
Assume I sample a signal according to the Nyquist criterion. Then I perform a simple summation / integral over a linearly interpolated signal. Is this equivalent to the integral over the continuous signal, i.e.?

$$\sum_{n=-\infty}^{\infty}x(nT) T = \int_{-\infty}^{\infty}x(t)\text{d}t$$

I hope my representation is correct for what I mean. My question relates to whether a reconstruction, e.g. Whittaker–Shannon interpolation, is necessary   in practical application. I would also appreciate any literature pointers for discussions on how this might affect practical cases with finite time and imperfect filters.


texts:
The expression for Fourier transform is given as
$$ F(\omega) = \int\limits_{-\infty}^{+\infty} f(t) \ e^{-j \omega t} \ dt  \tag{1}$$
Now, let one function be $f(x)$ and other be $e^{j\omega t}$ then their convolution at $p(0)$ is given as
$$ p(0) = \int\limits_{-\infty}^{+\infty} f(t) \ e^{-j \omega t} \ dt  \tag{2}$$ same as in $(1)$.
Also, if we split the exponential term into sinusoidal,the expression $(1)$ becomes
$$ F(\omega) = \int\limits_{-\infty}^{+\infty} (f(t) \cos{\omega t} -j f(t)\sin{\omega t})\ dt  \tag{3}$$
Let $f(t)$ and $\cos{\omega t}$ be 2 functions. So their correlation  at value
$q(0)$  is given as
$$ q(0) = \int\limits_{-\infty}^{+\infty} (f(t) \cos{\omega t}) \ dt  \tag{4}$$ same as in $(3)$. (we can do with f(t) and sine function also.)
So I want to know

Whether Fourier transform formula be considered as Convolution or Correlation?

Also, what is meaning of $p(0)$ and $q(0)$ terms ?


texts: temporal smoothing: FWHM of Gaussian kernel vs. window length of moving average

texts:
Suppose in one case I convolve Gaussian kernel with FWHM=10 samples. I would like to compare the result with moving average. My question: should I take the moving average window also 10 samples? In other words, is the size of FWHM equivalent to window of moving average?

Thanks!


texts: Re sizing and thresholding small font scanned Arabic words

texts:
I am working on an arabic OCR for printed scanned documents ... Some of the scanned documents may written with font size of height 8 which is quite small ... I want to resize the height to 60 pixels but some artifacts may occur due to nature of arabic characters .. some characters may overlap .I had used local thresholding methods after resizing but the results are still not acceptable ... any ideas ?

This is an example image:

This is the same example after resizing and applying local adaptive thresholding using 50 as a window size:

As you can see there are some discontinuities in some characters like this one :

Is there is any method to resize the image while preserving textual shape ?

My approach to fix characters discontinuities :

Threshold original image before resizing using local adaptive threshold using window size of 16 (this will solve characters discontinuities but the holes in the characters are filled) call it smallbw.
Resize smallbw using imresize(smallbw, [nh nw], 'nearest') and fill the holes in characters using imfill
Resize the original image to height 60 pixels using imresize(originalIm, [nh nw], 'nearest') call it largebw
Fill holes in largebw using imfill and call it bwfill
Extract holes from largebw by bwholes = bwfill - largebw
Finally, subtract bwholes from smallbw to get this

you can see here that the discontinuity found in the character @Image 3 has been solved ... but there is another problem raised here some characters may overlap as shown here

This is the best results I could achieve so far ... are there any other ideas that may solve these issues ??
and if you think this issue has no solution how would I resolve it rather than using resizing? what about using 12 font size text instead of 8?

Useful links :
Local adaptive threshold method used

Operating system : windows 7

Programming environment : Matlab 2013a - Image Processing toolbox


texts: Spectrogram from Complex Morlet wavelet ( recreating article result from data)

texts:
I am trying to recreate the following results:

...from the following data: Dropbox Link (FS40.dat)

The time is in milliseconds (frequency is 40000 Hz)

The article state that they used Complex Morlet wavelet to create the spectrogram:
"Power estimates from the averaged LFPs were calculated from
time–frequency spectrograms of the data from 1–88 Hz by convolving the signals with a complex Morlet wavelet of the form

$$
w(t,f_{0})=Ae^{(\frac{-t^2}{2\omega^2})}e^{(2i{\pi}f_{0}t)}
$$

for each frequency of interest $f_{0}$, where $\omega=\frac{m}{2{\pi}f_{0}}$, and $i$ is the imaginary unit. The normalization factor was $A=\frac{1}{\omega(2{\pi})^0.5}$, and the constant $m$ defining the compromise between time and frequency resolution was 7.

I only managed to get some "good" results using spectrogram function in matlab.
But I dont have much idea of how to use the morlet complex wavelet.
I got bad result when trying to use cwt with 'morl' window

Thank You.

P.S.

I'm trying to recreate this article:

Computational modeling of distinct neocortical oscillations driven by cell-type selective optogenetic drive: separable resonant circuits controlled by low-threshold spiking and fast-spiking interneurons.


texts: Classifying 2D shapes by the smoothness of their boundaries

texts: I'm doing image analysis, and I want to classify smooth objects (has smooth boundaries) from non-smooth objects (has zigzag-like boundaries). Which feature should be fed into ML framework? What are some popular techniques for feature extraction of shapes?


texts: When plotting the cosines of the phases of some sensors, I got an interesting plot

texts:
I have a laser that I modulate at some frequency, typically 1250 Hz. I have a sensor that tracks what happens inside a reactor when it's illuminated. I am doing Fourier analysis on both the laser modulation signal and the sensor, and extracting the components at 1250 Hz, forming a time series of phases.

I was playing around with these phase relationships in Excel and I came up with this plot. Laser is the x-axis, sensor is the y-axis. I'm trying to think back to college because I swear I've seen something like this before, but would appreciate any commentary!


texts:
Can we change the scale of image by using FFT? I mean, how should i do process on frequency domain of image to upscale or downscale the orginal image?

The other question of mine is that how can changing scale of image affect on fft of this image?


texts: Extended kalman filter for linear system

texts: I studied about Kalman filter(KF) for about a week. And I understand that Kalman filter(KF) is suitable for linear system and extended Kalman filter(EKF) can be used for nonlinear system. However, I wonder if extended Kalman filter(EKF) can be used for linear system. I understand that Kalman filter(KF) would be okay for the linear system. But I wonder if there's logical error of using the extended Kalman filter(EKF) for linear system.


texts: How to use low pass filter in Matlab for discrete data?

texts: As see in the picture, there are high frequency noise in the signal. How to design a low-pass filter to eliminate the noise? How to apply this filter in Matlab?


texts:
I have an accelerometer that I am reading the x, y, z values for m/s^2. I've written up my code in Python as follows.

while True:
  x, y, z = sensor.accelerometer
  with open('data.txt', 'a+') as f:
    f.write('{},{},{}\n'.format(x, y, z))
    f.flush()
    f.close()
  time.sleep(0.1)

The time.sleep(0.1) says to sleep for 1/10 of a second, so theoretically, in 1 second, I should have 10 samples. Thus, I would assume my sampling rate is 10 samples per second. However, what I have noticed is that I do not get 10 samples per second; there is overhead with the file writing (in actuality, I am doing more stuff than just reading and writing the accelerometer data). In fact, in one second, I get a distribution over the number of samples per second.

Now that I have recorded the data, I need to analyze these signals, but with a lot of the signal processing examples, I see that the sampling frequency is a required parameter.

My question is, what should I set as my sampling frequency?

Should I put in 10? If so, empirically, I am not consistently (if at all) getting 10 samples per second.
If I do not use 10, I have a distribution over the number of samples per second, should I just use the mean? If I use the mean number of samples per second as the sampling frequency, do I need to remove samples within a second interval?
What do I do about a second of interval with less than the mean number of samples?

Any help into these questions is appreciated.


texts: What is the effect on Hilbert transform of signal after multiplication by sinusoid?

texts:
Question :

Find Hilbert transform of $[u(t-a)-u(t-b)]\cos2\pi f_{0}t\\\\$ such that $\\0<a<b$

my attempt:

we know Hilbert transform of

$[u(t-a)-u(t-b)]\xrightarrow{\mathcal H} \dfrac{1}{\pi}\ln\left|\dfrac{t-a}{t-b}\right|$

but after multiplying left side by sinusoid $\cos 2\pi f_{0}t$  how the right hand side will alter??. Can we expect good results ( like shifting in spectrum by $\pm f_{0}$ which takes place in Fourier transform)


texts: How to decompose a mixed signal into its original signals?

texts:
I have generated a number of sinusoidal signals and added them together using the following Matlab code snippet:

%% Time Properties
Datapoints = 512;               % Sampling frequency
T = 1/Datapoints;               % Sampling period
Length = 512;                   % Length of signal
t = (0:Length-1)*T;             % Time vector

%% Signal Properties
signal=0;
for amp=1:10
    for freq=2:5:100
        for ph=0:pi/36:pi
            signal=signal+amp*sin(2*pi*freq*t+ph);
        end
    end
end

After that I get back the original frequencies using a fast Fourier transform:

NFFT=Datapoints;
fs=Datapoints/Length;
f = fs/2*linspace(0,1,NFFT/2+1);
fftsignal=fft(signal);
fax_bins=[0:NFFT-1];
[pks,pos]=findpeaks(abs(fftsignal.^2),fax_bins);

Now pos will give me the constituent frequencies.

In addition to that I need to know the corresponding phases and amplitudes to reconstruct the original signals.

How can I get the original phase and amplitude information?

Edit 1 But pks is not giving me the constituent amplitudes.

Edit 2 For debugging I used the following code snippet:

Fs = 20000;
t = 0:1/Fs:0.01;
fc1=200;
fc2=300;
x1 = sin(pi*fc1*t);
x2 = sin(pi*fc2*t);
x=x1+x2;
x=x';
xFFT = abs(fft(x))/length(x);
[amp,freq]=findpeaks(xFFT);

I get amp=0.6252 and freq=2


texts:
I have plot some of the IIR filter's impulse responses and some of them are exponentially increasing, increasing, decreasing etc but some IIR's impulse response approaches 0. If for example I have an IIR filter defined by its lccde:

$$ 6 y[n] + y[n-1] - 2y[n-2] = x[n] - 2x[n-1] $$

or

$$ y[n] = \frac{1}{6}x[n] - \frac{1}{3}x[n-1] - \frac{1}{6}y[n-1] + \frac{1}{3}y[n-2] $$

when I plot its impulse response in matlab:

[h,n] = impz([1 -2],[6 1 -2]);
stem(n,h);

its impulse response would approach 0, and this zero value will bestarting at n = 23. Does this mean that the IIR is an FIR since I can cut off the impulse response at n = 23?


texts: Example of radial frequency functions that give scale shiftability

texts:
Simoncelli's "Shiftable multi-scale transforms" paper talks about shiftability of the parameters of a wavelet transform, in particular translation, rotation and scale. The first two are used in his steerable pyramid decomposition.

Can anyone point to an example of a filterbank used for the latter - scale shiftability? I have used sin and cos (Fourier series) as radial basis functions, but I wonder if there is a better set that have been optimised for some property, e.g. space-frequency localisation.


texts: How to get the MLSEEqualizer work with IFFT and FFT blocks in the communication system?

texts:
Matlab MLSEE equalizer works fine with the QPSK modulated signal when the signal is propagated
over the multipath channel. However, if there are IFFT and FFT blocks used at TX and RX sides, respectively,
the receiver in the modified code has worse BER results than the original code as shown in the following MATLAB codes.

I would like to ask for your help to investigate the failure in the modified code below.

% Original Code: Works fine and gives perfect BER results.
%   Equalize a QPSK signal transmitted through a dispersive channel
%   using MLSE

qpskMod = comm.QPSKModulator(0,'SymbolMapping','Binary');
qpskDemod = comm.QPSKDemodulator(0,'SymbolMapping','Binary');
% Channel coefficients
chCoeffs = [.986; .845; .237; .12345+.31i];
N=512; % Modulated signal length
mleq = comm.MLSEEqualizer('TracebackDepth',10,...
                  'Channel',chCoeffs, 'Constellation',[1 1i -1 -1i]);
% Create an error rate calculator
ber = comm.ErrorRate;
for n = 1:10
  data= randi([0 3],N,1);
  modSignal = qpskMod(data);
  % Introduce channel distortion.
  chanOutput = filter(chCoeffs,1,modSignal);
  % Equalize the channel output and demodulate
  eqSignal = mleq(chanOutput);
  demodData = qpskDemod(eqSignal);
  % Compute BER
  a = ber(data, demodData);
b=a(1)
end

% Modified Code: IFFT and FFT blocks are used but the system shows worse BER results.
%   Equalize a QPSK signal transmitted through a dispersive channel using MLSE

qpskMod = comm.QPSKModulator(0,'SymbolMapping','Binary');
qpskDemod = comm.QPSKDemodulator(0,'SymbolMapping','Binary');
% Channel coefficients
chCoeffs = [.986; .845; .237; .12345+.31i];
N=512; % Modulated signal length
mleq = comm.MLSEEqualizer('TracebackDepth',10,...
                  'Channel',chCoeffs, 'Constellation',[1 1i -1 -1i]);
% Create an error rate calculator
ber = comm.ErrorRate;
for n = 1:10
  data= randi([0 3],N,1);

 % Modulate the data and convert it to time domain.
  modSignalx = qpskMod(data);
  modSignal=sqrt(N)*ifft(modSignalx);

  % Introduce channel distortion
  chanOutput = filter(chCoeffs,1,modSignal);

  % Equalize the channel output and demodulate
  eqSignalx = mleq(chanOutput);
  eqSignal=(1/sqrt(N))*fft(eqSignalx);

  demodData = qpskDemod(eqSignal);
  % Compute BER
  a = ber(data, demodData);
b(n)=a(1);
end

Many thanks in advance!

@Zeyad_Zeyad I have included the whole code here. It gives high error rate even the SNR is set to 30.

clear all;
close all;
qpskMod = comm.QPSKModulator(0,'SymbolMapping','Binary');
qpskDemod = comm.QPSKDemodulator(0,'SymbolMapping','Binary');
% Channel coefficients
chCoeffs = [.986; .845; .237; .12345+.31i];
N=512; % Modulated signal length
mleq = comm.MLSEEqualizer('TracebackDepth',10,...
                  'Channel',chCoeffs, 'Constellation',[1 1i -1 -1i]);
snr=30;
numiteration=10;
for n = 1:numiteration
  data= randi([0 3],N,1);
dataIn(:,n)=data;

% Modulate the data and convert it to time domain.
modSignalx = step(qpskMod, data);
modSignal=sqrt(N)*ifft(modSignalx);

% Introduce channel distortion
chanOutputx = filter(chCoeffs,1,modSignal);

    % Add awgn noise
 chanOutput=awgn(chanOutputx,snr,'measured');

    % Equalize the channel output and demodulate
    eqSignalx = step(mleq, chanOutput);
    eqSignal=(1/sqrt(N))*fft(eqSignalx);

    demodData = step(qpskDemod, eqSignal);
demodDataOut(:,n)=demodData;
    % Compute BER
  end
a = biterr(dataIn(:), demodDataOut(:));
b=a/(N*numiteration)`


texts:
The problem is the non- negative matrix factorization of a matrix. Let me explain my problem I have an original matrix $A=\begin{bmatrix}
0.248437 &0.25198098 & 0.25396825 & 0.25077881\\
0.253125& 0.25198098& 0.24761904 &0.24766355\\
0.248437& 0.24722662& 0.24920634& 0.24766355\\
0.25& 0.24881141& 0.24920634& 0.25389408
\end{bmatrix}$ and I have a matrix $B=\begin{bmatrix}
0.38331781& 0.385229& 0.384043 &0.381118\\
0.37776086& 0.377539& 0.381224 &0.381850\\
0.29993458& 0.300691& 0.298606 &0.299959\\
0.34108473& 0.337021& 0.335967 &0.337317
\end{bmatrix}$. What I know is that the matrix $B$ is some scaled version of $A$ i.e $PD$ where $P$ is a permutation matrix and $D$ is a diagonal matrix. Can i somehow get an approximation of $A$ from $B$ if I retain some properties of $A$ like the maximum of each row and minimum of each row to get an idea of how much scaling has been done, similarly if I keep a record of each column sum , and its maximum and minimum value. My main question is can I estimate my $A$ from $B$ by keeping some parameters with me beforehand. Can somebody help or suggest some text for this.


texts:
Is there any way one can prove that $\text{rank}\big(F+F^{T}\big)=3$? where F is a fundamental matrix. If you have any idea, do let me know.
I am trying to solve this problem.

I am just trying to understand the first part.


texts:
Let's say I have the very same signal stored in 2 .wav files. And lets say I process each file using a function, and it is established that one file is 20 phons larger than the other. Lets ignore how that function calculates this.

Does this mean that one file is approximately 4 times louder than the other?

What formula can I use to relate the difference in phons for two signal to the difference in their loudness?

-10 phons = 0.5
0 phons = 1
10 phons = 2


texts:
I'm having trouble showing the following relation:

The autocorrelation of a LTI system with impulse response $h[n]$, input $x[n]$ and output $y[n]=h[n]*x[n]$ is given as:

$r_{yy} = r_{hh}*r_{xx}$

My approach was the following:

$r_{yy} = E(y[n]y[n+m]) = E((h[n]*x[n])\cdot (h[n+m]*x[n+m])) = E((h[n]h[n+m])*(x[n]x[n+m])*(h[n]x[n+m])*(x[n]*h[n+m])) $

but I don't know how to simplify this further...


texts:
During Cross-correlation process:

xcorr=ifft(fft(a)fft(b)∗)

we have frequency offset and code phase offset.

How is it used later in the process of finding a user location?


texts: Spectral Correlation Function estimation using FAM method -- implementation help

texts:
I am trying to estimate the Spectral Correlation Function by the FFT Accumulation Method (FAM) in Matlab.  The way I understand the FAM method is as this link (see page 6 & 7 for FAM tutorial).

Here is the image again summarizing the process:

I'll summarize how I am understanding this method below:

Acquire a data chunk, of size = N data points.
Have the user enter in the desired frequency resolution (d_f), and desired alpha resolution (d_alpha).
From d_f and d_alpha, calculate the following parameters.  N' = fs/d_f; fs is sampling frequency of input signal.  N' turns out to be the STFT FFT size.  L = N'/4, where L is the STFT skip (see diagram in link).  P = fs/(L*d_a), which is the smoothing FFT size .. and finally N = P*L, N is the total number of data points we are going to process.  N should be nextpow2 > the input size and zero padded if necessary.

I have the code to chunk up the data into segments of N' data points, perform windowing for each of those chunks, compute the FFT of each of those chunks, and perform the downshift to baseband.  I store the complex baseband into a matrix called XD.  The dimensions of XD are: N' rows X P cols, where each col represents the STFT for a chunk.  Then, I perform the following operations:

XD = conj(XD');  % now dim(XD) = P x N'.  we do this b/c FFT operates on columns in a matrix in matlab.

The steps I am confused about in the FAM method are the multiplication of X with its complex conjugate, performing the time smoothing w/ the P point FFT, and then outputting the SCF.

To perform the multiplication w/ the complex conjugate, I do the following:

% Np = N' (in code).
XM = zeros(P, Np^2);
for k=1:Np
    for l=1:Np
        XM(:,(k-1)*Np+l) = (XD(:,k).*conj(XD(:,1)));
    end
end

Essentially this operation is creating a matrix XM with dimensions of P x (N')^2.  After this, I do an FFT on XM

XM2 = fft(XM);

Now what do I do?  I have to essentially select certain values from XM2 as my SCF estimate, but am confused on how to go about doing this?  Each row in XM contains the spectral correlation between different frequency bins.  So for example, column 1 of XM contains the spectral correlation of bin1 with bin1, col 2 contains spectral correlation of bin1 with bin2, col N' contains the spectral correlation between bin1 and bin-N', col N'+1 contains the spectral correlation between bin2 and bin1, etc etc...  That is the reason that XM has dimensions of P x (N')^2.

I'm confused as to which of these points I need to select in order to create my SCF estimate.  I essentially need to figure out which elements of XM correspond to a certain alpha and freq value, but I don't know how to do that.

Thanks in advance, and happy new years!


texts: Erosion compotition rule $\left ( f \ominus g \right ) \ominus h = f \ominus \left (g \oplus h \right ) $

texts:
I want to prove that erosion of a signal by a function $g$ followed by erosion with another function $h$ is equivalent to erosion of the signal by the dilation of the two functions (Erosion property):

$$\left ( f \ominus g  \right ) \ominus h = f \ominus \left (g \oplus h  \right ) $$

My approach:

By definition of erosion we get that: $$\left ( f \ominus g\right )(x) = ⋀_{y} f(x+y)-g(y)$$

So, we get that $$\left ( f \ominus g  \right ) \ominus h = ⋀_{y}\left ( ⋀_{z} f(x+y+z)-g(y)-h(z)\right )$$

Is my approach correct so far? And if yes, how do I go on?


texts: Deconvolution of shifted gaussian function in the frequency range

texts:
I have a signal defined as
$$A(t)\cdot\exp\left(-i\omega_0t\right)$$
with $A$ the envelope function and $\omega_0$ the carrier frequency. I would like to transfer this signal into the fourier space and back, but need a large amount of points in time to accurately resolve the carrier frequency (the larger the time window, the more points are necessary). To circumvent that issue I intend to use a different strategy by using
$$F(A(t)\cdot\exp\left(-i\omega_0t\right))=\widehat{A}(\omega)\star\delta\left(\omega - \omega_0\right)$$
After the amount of points necessary to resolve $A$ accurately is significantly smaller than to resolve $\omega_0$ that significantly speeds up my calculation. Nevertheless, to apply the inverse Fourier transformation with the original amount of points used for the transform from $A(t)\rightarrow\widehat{A}(\omega)$ I have to "deconvolve" this function again by convolving it with
$$\delta(\omega + \omega_0)$$
but would require having negative frequencies. Is that even possible? If not, are there other ways to "shift" the function back, such that I can use fewer points for my FFT?


texts: Forward error correction with limited possible errors

texts:
I am using a Reed-Solomon FEC implementation in GF(1024) to correct T=15 errors in a given codeword.

In the particular implementation though, I know for a fact that for each transmitted symbols, the corresponding received (potentially error'ed) symbol will always be one of a given set of symbols. The set depends on the value of the original symbol, and can contain anywhere from 32 to 243 members.

So for example if my original symbol is A0='0000000000', then the received symbol will always belong to a known set S(A0) which obviously contains A0 and some other symbols, never to exceed 243 symbols total per set. If A!=B then S(A)!=S(B) (but the intersection of S(A) and S(B) in general is non-empty).

Can you think of some optimization that I can use in order to reduce the parity of the code, given that I know that there is a limited number of possible received symbols for each transmitted symbol?

Thanks,
-Dimitri


texts: how to create a harmonic mask from fundamental frequency?

texts:
Using an algorithm that gives the fundamental frequency across time, example https://www.mathworks.com/help/audio/ref/pitch.html (or some other one?)

how would I then use that, to create a harmonic mask of nHarmonics based on that fundamental frequency? So the mask will contain the f0 + nHarmonics harmonics.

The point of this mask will be to then use it as a filter for the audio mixture.

I want to do this using MATLAB.

edit: I forgot to mention that the f0 is not a constant value, it is a ranged f0 track from fmin to fmax, for example dominant f0 from 100Hz to 800Hz. So the f0 track contains different values, not just one value. I can understand a comb filter working for one value, but this is a bit more complex, as the f0 track has multiple different values, so now each one of those would have to have nHarmonics multiples above the f0 for me to create this harmonic mask.


texts: Is it correct to subtract a low-pass filtered signal from the original signal and use the result as a "high-pass"?

texts:
I am having a hard time trying to find documentation to implement band-pass or high-pass filter with python/scipy/numpy.

I can easily create and apply a low-pass filter, though, so I ask:

Would it be conceptually correct to low-pass-filter a signal, then subtract the result from the original signal, in order to get just the high-frequencies?

Also, if anyone has a simple example of a naive bandpass filter in Python (preferrably using the numpy and scipy libraries), I'd be very thankful.

What I look for is something like:

filtered_signal = band_pass(original_signal, rate, low=20, high=500)

Thanks for any help!

EDIT: with scipy, I'm using this as low-pass, with good results:

import numpy, scipy.signal

def firfilt(interval, freq, sampling_rate):
    nfreq = freq/(0.5*sampling_rate)
    taps =  sampling_rate + 1
    a = 1
    b = scipy.signal.firwin(taps, cutoff=nfreq)
    firstpass = scipy.signal.lfilter(b, a, interval)
    ## second pass to compensate phase delay
    secondpass = scipy.signal.lfilter(b, a, firstpass[::-1])[::-1]
    return secondpass


texts:
I'm a newbie to Signal Processing - my apologies if this question is too obvious (I'm a financial trader trying to use DSP techniques).

For a linear filter:

$y[n] = (1-p) x[n]+p y[n-1]$

we can the write Transfer Function as

$H(z) = \frac{1-p}{1-pz^{-1}}$

But instead, suppose I have a non-linear filter:

$y[n] = y[n-1]+(1-p)\frac{1}{2}(x[n]e^{-2 y[n-1]}-1)$

How can I write a Transfer Function ?

For example, if I simply replace $y[n-1]$ by $y[n]z^{-1}$ and rewrite I get:

$H(z) = \frac{(1-p)\frac{1}{2}(x[n]e^{-2 z^{-1} y[n]}-1)}{1-z^{-1}}$

which is a function of $x[n]$ and $y[n]$ in addition to being a function of $z$?

Can this be written merely as a function of $z$?


texts: Fourier Transform: $\omega$ vs $f$ as frequency variable

texts:
I try to understand how the Fourier transform changes when I try to compute $X(\omega)$ or $X(f)$. Can someone work me through the maths please for two examples, $x(t) = \exp(j\omega_0 t)$ and $x(t) = \cos(\omega_0t)$.
I know what their Fourier transform is in terms of $\omega$ but when I try to compute $X(f)$ from first principles, I am stuck.


texts: I cannot find any reliable or reasonable differences between the two even after searching for it! Please tell me the possible differences with reasons and reference links!


texts:
I have a sensor that applies a frequency dependant phase alteration to my signals.  I'm trying to find a way to correct these phase responses using an FIR filter.  Essentially, I'm looking for methods to apply an all-pass filter with arbitrary frequency response to correct this.
For magnitude calibrations of my sensor, I can achieve this by drawing out my calibration spectrum, applying an ifft to it, then applying a hamming window.  Then I simply convolve the resulting impulse response with my signal.  This is essentially MATLAB's fir2().  Is there an analogous procedure I can do to correct arbitrary phase responses in my signals?